Supersolidity of cold-atom Bose-Fermi mixtures in optical lattices 
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We investigate a cold atomic mixture of spinless bosons and fermions in two-dimensional optical 
lattices. In the presence of a nested Fermi surface, the bosons may develop a fascinating supersolid 
behavior characterized by a finite superfluid density as well as a spatial density wave order. Focusing 
on the triangular lattice geometry and combining a general Landau-Ginzburg- Wilson approach with 
microscopically derived mean-field theory, we find an exotic supersolid phase at a fermionic band- 
filling of n/ = 3/4 with a Kagome-type crystalline order. We also address the case of anisotropic 
hopping amplitudes, and show that striped supersolid phases emerge on the square and triangular 
lattices. For weak interactions, the supersolid competes with phase separation. For strong intra- 
and inter-species interactions, with the total number of fermions and bosons corresponding to one 
particle per site, the bosons form an alternating Mott insulator ground state. Finally, for a mixture 
of «^Rb«K and ^^Na^Li, we show that supersolidity can be observed in the range of accessible 
temperatures in the square lattice geometry. 

PACS numbers: 03.75.Lm; 67.80.kb; 67.85. Pq 



I. INTRODUCTION 



One of the most intriguing predictions in the theory 
of quantum mechanics is the possibility of supersolidity 
- superfluid behavior in a rigid crystal (solid). A super- 
solid phase involves two unrelated broken symmetries - 
global U{\) phase invariance breaking (superfluidity) and 
translational invariance breaking (density wave) P, |2l ■ In 
the context of a lattice system, the discrete translational 
symmetry will be broken (forming a superlattice struc- 
ture). Over the years, much work has been devoted to its 
experimental realization as well as its theoretical under- 
standing 0, H, H, Q . While experimental efforts have so 
far concentrated on solid "^He, there is a variety of the- 
oretical models that exhibit supersolidity, most notably 
interacting lattice models such as the Bose-Hubbard 
or the Bose-Fermi Hubbard model. The possibility 

to simulate these models u sing ultracold atoms in opti- 
cal lattices [HI [H [H, Q [il m, [13, H [la, 20] offers 
a fascinating alternative experimental route to superso- 
lidity [m, [H, [H [^, [Hi. In particular, time of flight 
experiments can probe superfluid condensation and den- 
sity wave order simultaneously [1^ . 

Single-component Hubbard models require a nearest- 
neighbor interaction, at least, to stabilize a supersolid 
phase Izl, [IS m [13, [U, [13, [H, [H, H; a typical ex- 
ample being the frustrated bosonic U-V model - where 
U denotes the on-site Hubbard interaction and V the 
interaction between nearest neighbors - on the triangu- 
lar lattice. Another concrete example is the dipolar bo- 
son lattice model [sgI W\ . On the other hand, it has 
been shown that in two-component mixtures already pure 
on-site interactions are sufflcient to induce a supersolid 
phase [2l|, [H, [H, [i^]. Mixtures with different species, 
in general, allow to realize new states of matter in a 
variety of settings El,!!!, [H, [11, [H, [H, S [H, [H, 
[sol . [Sll . [5^ ISSl . kA ISSI . l56j . due to effective interactions 
that one species induces on another. Hereafter, we fo- 



cus on Bose-Fermi mixtures in which the fermions in- 
duce an effective interaction between the bosons, and 
vice versa. We study in detail the emergence of su- 
persolidity in the Bose-Fermi mixtures; other aspects of 
Bose-Fermi mixtures have been addressed in the litera- 
ture [13, [11, [li, [10, [Ml, [11 [11 [el [Mils [11 [Mill [23 ■ 

In Bose-Fermi mixtures, one important possible mech- 
anism to achieve boson supersolidity relies on the exis- 
tence of a nested Fermi surface. With nesting, fermions 
tend to exhibit a density wave at the nesting wavevec- 
tor(s), and this generates the same ordering tendency 
on the bosons, through the boson-fermion interaction. 
Alternatively, the fermions induce interactions between 
bosons, and the superfluid-to-supersolid transition can 
also be understood as a condensation of rotons, oc- 
curing when the roton gap in the superfluid excita- 
tion spectrum vanishes upon increasing the interaction 
strength [2^, H^, |7l| . We will employ two distinct mean- 
field calculations which follow both points of view. 

Here we focus on two-dimensional lattices with both 
spatially isotropic and anisotropic hopping amplitudes; 
the anisotropic cases will allow us to investigate the 
(quasi-)one-dimensional limit of supersolidity. In one 
dimension, where fermions are formally equivalent to 
hard-core bosons, supersoHdity was recently predicted to 
occur in a strongly interacting two-component bosonic 
mixture (40| . also as a non-equilibrium state [23| . In 
three dimensions, our analysis predicts that the super- 
solid appears at lower temperatures compared to the two- 
dimensional case, because the presence of van Hove sin- 
gularities strongly enhance the tendency for density wave 
formation only in lower than three dimensions. 

We begin our analysis with a mixture of spin-polarized 
bosons and fermions in a two-dimensional triangular op- 
tical lattice [T^I, which exhibits nesting at a particular 
fermionic band filling of n/ = 3/4, as shown in Fig. [1] 
A supersolid phase has already been predicted for the 
isotropic square lattice [H HI] . For a sufficiently deep 
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optical lattice only nearest-neighbor hopping survives, 
and the system is described by the ubiquitous single band 
Bose-Fermi Hubbard Hamiltonian [l^, , 

H = -J2^f,^Jfjf,+tt^^Ab,) (1) 

- {fifmj + fibUj) 

i 

+ '^'("^ ~ 1) + ^b/ "^i^^ > 

where f^{bl) is the fermionic (bosonic) creation opera- 
tor at site i, while nii = fjfi (rii = b\bi) denotes the 
fermionic (bosonic) number operator and is the 

chemical potential of the fermions (bosons). Hopping is 
restricted to neighboring sites (denoted by the summa- 
tion over with amplitudes tf,ij{tb^ij) for fermions 
(bosons) that in general depend on the direction of hop- 
ping. The on-site boson-boson and boson-fermion inter- 
action strengths are given by Ubb and Ubf, respectively. 
We do not include interactions between fermions, which 
due to the Pauli-principle only occur in the p-wave scat- 
tering channel, and are frozen out at ultracold tempera- 
tures. 

Before we set out to show how and in which parameter 
regime supersolidity arises from the above Hamiltonian, 
we give a brief outline of the article. 

The paper is organized as follows: after this introduc- 
tion, in Sec. [TTl we focus on the isotropic triangular lat- 
tice, and carry out an instability analysis of the system in 
the weak-interaction limit after the fermionic degrees of 
freedom were integrated out. For the particular fermionic 
filling rif ~ 3/4, the Fermi surface both shows nesting 
and contains van Hove singularities. This triggers two 
distinct low-temperature instabilities of the superfluid 
bosons: one towards phase separation and one towards 
supersolid formation. The density modulation in the su- 
persolid phase is characterized by the three nonequivalent 
nesting vectors Qi_2 3 of the Fermi surface (see Fig.[T]^b)), 
producing a Kagome-type crystalline order in real-space. 

In Sec. IHI[ we calculate the low-temperature mean- 
field phase diagram of the system, containing a Kagome- 
supersolid and a phase separated regime. We also de- 
termine the amplitude of the density wave modulations 
inside the supersolid. 

In Sec. IIII Al we use a general Landau-Ginzburg- 
Wilson mean- field theory to find a number of likely super- 
solid phases in a general phase diagram, assuming pos- 
sible condensation into the wavevector modes 0,Qi^2,3- 
We match the phenomenological Landau expansion pa- 
rameters to a microscopically derived mean-field expres- 
sion and obtain criteria where the supersolid and phase 
separated regions emerge in the phase diagram, taking 
us beyond the instability analysis. 

We then embark, in Sec. IIIIBl to calculate the density 
wave modulation in the supersolid using a different mean- 
field approach that treats the fermions exactly. We find 



that the fermionic spectrum acquires a gap in the super- 
solid phase, allowing the system to lower its energy. The 
density modulation in the supersolid is found to be rather 
weak, typically involving only 0.1% of all the bosons. The 
transition temperatures Tss are also small compared to 
Tp, typically we find Tss/^f — 0.01, where we define 
Tp = &tf as the Fermi temperature of the optical lattice 
(the Boltzmann constant fc^ — 1). 

To find larger transition temperatures, we turn to in- 
vestigate the case of anisotropic hopping amplitudes in 
Sec. IIVI first for the square (Sec. IIV Ap and then for the 
triangular lattice (Sec. IIV Bp . As a result of reduced 
symmetry, only one nesting vector occurs and the su- 
persolid phase exhibits a striped pattern in real-space. 
Since the relevant features of anisotropic hopping are al- 
ready captured by the square lattice geometry, we mainly 
focus on this case, where analytical results can be de- 
rived. Since nesting is fulfilled for a larger fraction of 
wavevectors, we find larger supersolid transition tem- 
peratures in our mean-field analysis. They are as high 
as Tss ^ tf {3tf/5) ~ Tf/4 {Tf/5) for the isotropic 
(anisotropic) square lattice. The supersolid density wave 
now involves up to 20% of all the bosons. 

In Sec.|Vl we predict the supersolid parameter regime 
for two experimental realizations of Bose-Fermi mixtures, 
®'^Rb'*°K and ^^Na^^Li. We show that while the transition 
temperatures for the triangular lattice geometry are be- 
yond current cooling limits, the supersolid phase on the 
square (isotropic and anisotropic) lattice, should be ac- 
cessible with current technology. We give optimal choices 
of experimental parameters in order to maximize Tgg. 
Anisotropic hopping offers the possibly crucial advantage 
to significantly weaken the tendency towards phase sepa- 
ration while retaining supersolid transition temperatures 
close to current experimental limits. 

We discuss in detail how the supersolid phase can be 
detected using time-of-flight measurements, and conclude 
that the detection becomes easier for a smaller ratio of 
bosonic to fermionic hopping amplitudes tb/tf, i.e. slow 
bosons. 

Finally, in Sec. IVIl we address the limit of strong in- 
teractions, where the system, at total unit filling (bosons 
and fermions combined), can be described by a quantum 
Heisenberg Hamiltonian with an additional gauge field 
arising from the celebrated Jordan- Wigner transforma- 
tion in two dimensions. We show that the bosons are lo- 
calized in an (alternating) Mott insulator phase, and the 
fermions feature a density wave with wavevectors equal 
to the nesting vectors. 

We summarize our results in Sec. IVIII and leave the 
details of a number of our calculations to the appendices. 



II. LOW TEMPERATURE INSTABILITIES 

In general, one expects a supersolid phase to occur 
for weak-interspecie interaction Ubf, since for larger in- 
teractions the mixture either phase separates or enters 
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FIG. 1: (Color online) (a): Triangular lattice in real space 
with our conventions of unit length Bravais lattice vectors 
{ai,a2} and hopping amplitudes {ti,t2}. 
(b) Reciprocal lattice vectors {Gi,G2} (dotted), first Bril- 
louin zone (thick hexagon) , and the Fermi surface at jj,f = 2tf 
{uf — 3/4) (thin inner hexagon). The Fermi surface exhibits 
three nonequivalent nesting vectors Qi , Q2, Q3 (solid arrows), 
which when folded back into the first Brillouin zone occur at: 
Qi = G2/2, Q2 = Gi/2, and Q3 = (Gi + G2)/2. They 
coincide with the critical points (corners) of the Fermi sur- 
face hexagon, which give rise to a van Hove singularity in the 
density of states at this filling (see Fig. 2(b)). 



Gi = 27r (^1, -^j and G2 = 2tt (^-1, . The real space 

lattice and the first Brillouin zone are shown in Fig. [TJ 

The Hamiltonian in momentum space then reads H = 
Hb + Hf + Hbf with 

qeBZ ki.ka 
qeBZ 

q.ki.kaGBZ 

(2) 

where the dispersion relation for the fermions (bosons) 
on the triangular lattice reads 

C/(fc)(q) = -Ai/(b)-2i/(6)iC0sgi-%(t)2C0Sy cos^— . 

(3) 

The hopping amplitudes t/(6)2} describe hopping 

along the direction ±(a2 — ai) and along ±ai 2, respec- 
tively. Most generally there are three different hopping 
parameters for the three directions on the triangular lat- 
tice, but we will only consider two of them being differ- 
ent. This already includes the interesting cases of weakly- 
coupled one-dimensional chains (i/(6)2 ^ ^/(b)i) ^^'^ 
transition to the square lattice { tf[b) 2 S> if(b)i)- Both 
scenarios will be discussed in Sec. IIVI Until then, we as- 
sume isotropic hopping amplitudes = tf(b)2 = ^/(6)- 
The imaginary time partition function of the system 
is quadratic in the fermionic degrees of freedom, which 
can therefore be integrated out exactly. 



a Mott-insulating state (large Uhb and Uhf) [3l]. We 
will separately address the strongly interacting regime 
in Sec. IVII For now, we focus on sufficiently weak boson- 
fermion interactions Ubf (we will specify the exact con- 
dition below). 



A. Definitions for the triangular lattice 

We begin by integrating out the fermionic degrees 
of freedom in an imaginary time functional integral ap- 
proach. Beforehand, it is convenient to write the Hamil- 
tonian of Eq. ([T]) in the Bloch state basis, where the an- 
nihilation operators read fj = SqsBZ /q^*'* ■ and 

— — ^qSBZ^q'^ 



' ~ VNl ^qeBZ ■ 
Here, the summation is over 



the first Brillouin zone (BZ) of the triangular lattice, Xj 
is the real-space vector to lattice site j and Nl is the 
number of unit cells. Our convention of unit length Bra- 
vais lattice vectors is ai = ^i, and a2 = (^^hi i 

i.e., |ai^2| — 1; for simplicity, the lattice constant is fixed 
to one. The reciprocal lattice vectors are then given by 



B. Effective bosonic theory 

Integrating out the fermions, yields formally 

The variable q contains a momentum and an imaginary 
time component q — (t, q), and the bare action derives 
from the respective parts of the Hamiltonian: 



Sb 



S 



q 

Jo 



b;^b, + Hb{b;,b,) 



d_ 



+ 0(q))/< 

q 



(5) 



q,ki,k2 
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Here [3=1 /T where T is the temperature and impHcitly 
all the fields have the same imaginary time component. 
Integrating out the fermions, we get a determinant de- 
pending on the boson density 



rrcff c 
Of, — Dh — 



5,-Trln (^g-\q)S{k - q) + ^ 



(6) 

Next we expand to second order in U^f. To first order 
in Ubf, the fermions simply produce a (trivial) shift of 
the bosonic chemical potential fJ-b fJ-b ~ UbfUf that 
depends on the fermionic filling Uf = Nj /Nj^ , where Nf 
is the total number of fermions. Trading the integration 
over imaginary time with a summation over the bosonic 
Matsubara frequency domain, defined by 



1 

7^ 



oo 

E ' 

m— — oo 



(7) 



with u!,n = 27rm//3, and simplifying a bit, the effective 
bosonic action up to second order in Ubf takes the form 

Sf^ E {Hu;„,+UcL)+Ubfnf]b;b, 

I3=(iw™,q) 

+ 2iO E ^^bb + UbfXiT,q)]bl^-gbl^+qbk^bk,Y 



ki,k2 



(8) 



The second order term in Ubf, involves the fermionic po- 
larization Lindhard function, 

which depends on temperature T via the Fermi func- 
tion /(^/) = [1 + exp(^//r)]~i. The induced interac- 
tion is attractive in momentum space independently of 
the sign of Ubf, as x(q) < for all q. In real-space, it 
is long-range and oscillatory in sign with an interesting 
(Kagome-lattice-type) structure due to the non-trivial 
wavevector dependence (see Appendix [B]) . 

Higher order terms can be neglected when MgUbf ^ 1, 
where Mq ^ tj^ is an estimate of the Lindhard function 
(away from its singularities; see Fig. (^^b)). 

We will analyze the behavior of this function in some 
detail in the following and in Appendix |^ since it pro- 
vides a basic understanding of the mechanism of super- 
solid formation. In particular, its static part {iLUm — > 0) 
diverges logarithmically for low temperatures at q = if 
the Fermi surface contains van Hove singularities. Fur- 
thermore, it diverges at special wavevectors q = aQi if 
the Fermi surface is nested with nesting vectors Q.;, and 
<C a < 1. If both features are present, as it is the 
case for a fermionic chemical potential ol fj,f = 2tf or 
fermionic band filling of = 3/4 (see Fig. [IJb) and 
Fig. (Ha)), the divergence at and close to the nesting 
vectors gets enhanced to x ~ where c is 




(b) 




6 -4 -2 2 3 

e (units of t,) 



FIG. 2: (Color online) (a): Fermi surfaces for f^f/tf = 2 
(dashed) and /x/A/ = {-3.75, -0.375, 0.875, 2.5} (solid). The 
thick hexagon denotes the first Brillouin zone, 
(b) Fermionic density of states g{e), which diverges logarith- 
mically at energy e = 2tf (van Hove singularity). We take 
Mo = S/'iiT^tf as a measure of (7(e) away from the divergence. 



a numerical constant. These divergences provide two 
competing low temperature instabilities in the superfluid 
(bosonic) phase, one towards supersolid formation and 
one towards phase separation. 

First we note, that it is legitimate to only consider the 
static limit {icUm 0) of x, if the fermions are much 
faster than the bosons {tf ^ tb). Then, the fermionic re- 
sponse occurs on much faster timescales than the move- 
ment of the bosons, and one can safely neglect retarda- 
tion effects. More formally, the terms with nonzero Mat- 
subara frequencies {ioJm 7^ 0) only contribute subdomi- 
nantly to the divergences at q = 0, Qi,2,3. The opposite 
limit of "slow" fermions, where superfluid bosons induce 
an attractive interaction among the fermions leading to 
(exotic) superconducting phases, has been discussed in 
Refs.[50„i58,]. 

It was shown in Ref. (6l| that the static approximation 
always yields qualitatively correct results, because (— x) 
is positive definite. In general, it improves for smaller 
fermionic densities n/. Note that we take the spatially 
non-local nature of the induced interaction fully into ac- 
count. 

In the static limit, one can work with an effective 
Hamiltonian for the bosons which takes the form: 



ki,k2 



t/(T,q) 
2Nj. 



hi hi , hir^bh:, 

ki— q k2+q ^2 ki 



(10) 



with an interaction U{T, q) that is given by, 

UiT,q) = Ubb + U^fX{T,q), (11) 

where x{T,(l) = xiT,iuj = 0,q). As mentioned above, 
this perturbative form of the interaction is valid for 
MoUbf <C 1, where Mo = 3/(47r^t/) is a measure of the 
Lindhard function and the density of states away from 
its singularities (as shown in Fig.[2l^b)). 
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At temperatures well below the Kosterlitz-Thouless 
transition temperature Tktj the bosons form a quasi- 
condensate, i.e., a condensate with a fluctuating phase, 
described by a wavefunction of the form ^n(x)e*'^'^^^ 
One can diagonalize H^^ in the superfluid phase em- 
ploying the well-established Bogoliubov approximation 
&q = V^LiT-oS (q) + 6q where only terms up to quadratic 
order in the fluctuation operators b^^o are kept and uq 
represents the flnite superfluid density. This yields the 
spectrum of elementary excitations of the superfluid 



Eb{T, q) = ^a(q)' + "^UbM (Ubb + U^fXiT, q) 

(12) 

Here, we have assumed that all the bosons are condensed 
into the zero momentum state by equating uq = rib, 
which is valid to a good approximation for temperatures 
T <gC Tkt- In the following, we discuss the two instabili- 
ties of the superfluid occurring when Eb (q) vanishes. 



C. Phase separation 



For small wavevectors |q| ^ 1, the Bogoliubov spec- 



trum is linear ^^(q) ^3nbtb[Ubb + U^fX{T,0)]\q\, 
with a sound velocity that vanishes at: 



xiT,0) 



^Ubb/Uf 



bf ■ 



(13) 



At this point, the contact interaction U{T, 0) becomes 
attractive, which marks the transition to a phase sepa- 
rated regime, since a Bose condensate is thermodynami- 
cally unstable for U{T, 0) < [zl. 

For a regular density of states g{e) at the Fermi sur- 
face, one finds that x{T,0) — —g{0). However, due to 
stationary points (|Vq^/(q)| ~ 0) on the Fermi surface 
for a chemical potential of /i/ = 2tf (see Fig. [2ljb)), the 
density of states diverges at this filling like 



8tf 



g{e) ~ Mo In 



where Mq — 3/(47r^i/), resulting in 



xiT,0) = -Mo In 



SC'itf 



T 



(14) 



(15) 



with Ci = 2e'-^/n w 1.13, and C being the Euler- 
Mascheroni constant. Thus, for any nonzero coupling 
Ubf between the bosons and fermions, there is a tem- 
perature Tpg''*- at which the q = term of the effective 
interaction becomes attractive [U{T^g^- ,0) = 0]: 



r^'r=8Cit/exp(- 



1 



BF 



(16) 



with Xbf = MoU^j/Ubb <C 1 describing the ratio of 
induced attraction to intrinsic repulsion between the 
bosons. 



D. Supersolid formation 

The other low temperature instability of the super- 
fiuid phase occurs only in the presence of a nested Fermi 
surface. Nesting is defined as the existence of a nesting 
vector Q such that for a finite domain of wavevectors k, 
the energy fulfills the prerequisite 



0(k + Q) = -C/(k). 



(17) 



Close to the Fermi surface with ^/(k) w 0, the denomina- 
tor in the expression of xi^, Q) becomes very small (see 
Eq. At the same time, the numerator is nonzero, 

since Q links an occupied with an unoccupied state, and 



/K/(k)]-/K/(k + Q)] 



0(k)-e/(k+Q) «/(k) 



-— • (18) 



Thus, nesting leads to the divergence of x{T 0, Q). 

On the triangular lattice, at the particular band filling 
of n/ = 3/4, as shown in Fig. [2^ a), the Fermi surface 
exhibits three nonequivalent nesting vectors which map 
the different sides of the Fermi surface hexagon onto each 
other. They read Qi = (— tt, 7r/\/3), Q2 — (7r,7r/\/3), 
Q3 = (0, 27r/-\/3), and coincide with the location of the 
van Hove singularities. They fulfill — Qi = Qi + G,„, 
with G„i being a reciprocal lattice vector, as well as Qi -I- 
Q2 = Q3 (and cyclic permutations). Each of them maps 
two of the six van Hove points onto another van Hove 
point, which leads to a significant enhancement of the 
divergence of x(T, Qi). 

We can analytically estimate this divergence of the 
Lindhard function (Eq. (O) by approximating. 



X{T,Q,) 



de 



g{e) tanh(e/2T) 
~3 2^ 



(19) 



where we have used the nesting relation ^^(k -f Qi) = 
— ^/(k), that strictly holds only for states along a rectan- 
gular path C, that goes in the case of Q3 along {— Qi ^ 
Q2 ^ Qi ^ — Q2 ~Qi}- We have inserted a fac- 
tor of 1/3, because the nesting property is only fulfilled 
for one third of the states on the Fermi surface. In ad- 
dition, we have ignored the fact that nesting is not ful- 
filled for all k-states while replacing the sum over the 
k-states that satisfy the nesting condition by an integral 
over all states. Nevertheless, solving the integral gives 



nesting vectors i 



^MSL 
6 



In 



SCitf 



which holds at the three 



= 1,2,3. If we compare this with nu- 
merical results using for example Monte-Carlo integra- 
tion, we observe that the slope Mo/6 is in perfect agree- 
ment, but the energy scale in the logarithmic function 
needs to be slightly adjusted. More precisely, we can 
easily fit our numerical results to the function: 



x(r,QO«--j^ 



Mo 
6 



In 



SaCit 



T 



(20) 



and obtain the fit parameter a = 2.17. 
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If we plug this result into the Bogoliubov dispersion 
relation (see Eq. (|12p ). we find that £'f,(Qi) decreases 
as the temperature is lowered, and finally becomes zero 
[Ebi^T■M^)=0] when 



bb 



Ui_ 
^bf 



^biQ^ 



2nbUbb 



(21) 



which defines the supersolid transition temperature 
based on the instability criterion 



rpinst. 

-'ss 



SaCitf exp 



3 



BF 



(2 + rB; 



(22) 



Here, tb = Stb/ni,Ubb is the ratio of kinetic to interaction 
energy of the pure boson system. 

The transition temperature Tgg*''' becomes larger for 
smaller tb, favoring slower bosons or larger intrinsic re- 
pulsion. However, one has to consider that at strong 
coupling Tb <C 1, there is a competing superfluid to 
Mott-insulator transition at commensurate densities n^, 
which occurs on the two-dimensional triangular lattice 
at the critical ratio {Ubb/'tb)c — 26.5 [73|. For a typical 
bosonic filling oi rib — 5/4, the Mott insulator appears 
at Tb ~ 1/4. In addition, weak-coupling requires that 
Ab_f < Tb which sets an upper limit to the value of Ubf. 

In short, the instability analysis provides an intuitive 
physical view on why we expect a condensation of rotons, 
i.e., (bQ-) 7^ 0, in the presence of a nested Fermi surface. 



Eq (units of tt) (a) 
1.2 




FIG. 3: (Color online) Bogoliubov dispersion relation along 
q = (gi,(/i/\/3) for various temperatures T and fixed param- 
eters ni,,ti,, C/;,;,, f/i,/. (a): Roton gap closes slightly away from 
the nesting vector Q2 = {■k,tv/\/3) for nj, = 1.25, tt — O.lt/, 
Ubb ~ "itf, Ubf ~ 1.9568f/ at the temperature logjgT/f/ = 
— 1.7 (lowest curve). Other curves correspond to the temper- 
atures logj^pT/f/ — —1, —1.3, —1.6 (top to bottom of upper 
three curves), (b): Roton gap closes at the nesting vector Q2 
for Ub = 1.25, tb = O.ltf, Ubb = 1.15t/, Ubf = l.llTtf at the 
temperature logjgT/tf = —2.9 (lowest curve). Other curves 
correspond to logj^Q T/tf — —1, —2.4, —2.8 (top to bottom of 
upper three curves). 



mainly focus on the commensurate supersolid phase that 
emerges below Tl (see Fig. |31Ib)). 

In the next Sec. lIIIi we study the phase diagram using 
more general bosonic and fermionic mean-field theories. 



III. PHASE DIAGRAM AND PROPERTIES OF 
KAGOME SUPERSOLID 



E. Incommensurate density wave 

It turns out that, at finite temperatures, one must be 
more careful with the analysis of the Lindhard function. 
We show in detail in Appendix IA 31 that in an intermedi- 
ate temperature regime, where the thermal smearing of 
the Fermi edge is larger than the level spacing (~ tf/N^), 
the minima of the Lindhard function occur at wavevec- 
tors slightly different from Q^. As a result, the roton 
gap closes (slightly) away from the nesting vectors, at 
Ki = aQi with a < 1 (see Fig. which leads to the 
formation of a density wave that is incommensurate with 
the lattice structure at intermediate temperatures. 

On the other hand, for a finite lattice composed of 
Nl unit cells, the level spacing starts to play a role at 
temperatures Tl ^ tf /N^ , and the minimum of the Lind- 
hard function shifts to at temperatures below T^, i.e., 
a — » 1 for T <^Tl, where thermal effects can be ignored. 
In this sense, the incommensurate regime does not sur- 
vive for T < Tl. A more quantitative analysis is given 
in Appendix I A 41 where we find that, 

Tl ~ 2TThf/NL. (23) 

Choosing an experimentally relevant lattice size of Nl — 
60 [7^, one finds log^g {Thltf) = —2.3. In this paper, we 
will not address in detail the properties of the (interme- 
diate) incommensurate density wave regime and we will 



In the following Sec. lIII Al we derive a low temperature 
phase diagram of the system using a bosonic mean-field 
theory that goes beyond the instability analysis. 

We identify a novel, highly symmetric supersolid phase 
with a Kagome-type density modulation in real-space, 
and calculate supersolid transition temperatures. In 
Sec. nil Bl we further study this supersolid and deter- 
mine the amplitude of the density wave modulation us- 
ing a different mean-field theory that treats the fermions 
exactly. 

A. Phase diagram from bosonic mean-field theory 

Here, we employ a Landau-Ginzburg- Wilson mean- 
field theory to build the low temperature phase diagram 
of the system. We find a novel Kagome-supersolid phase 
that competes with phase separation, and derive transi- 
tion temperatures to both phases generalizing the insta- 
bility resuhs of Eqs. HH) and ((22)l . 

1. Construction of the free energy 

Based on the results of the instability analysis, we ex- 
pect phase transitions to occur at low temperatures. We 
therefore construct a general Landau-Ginzburg- Wilson 
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free energy functional for the bosons on the isotropic 
triangular lattice. The details of this procedure can be 
found in Appendix [Cl 

We assume that the bosons may have a number 
of Fourier components condensing at momenta q = 
0, Qi,Q2,Q3- The complex bosonic order parameters 
{"00,1,2,3} are defined as 



(bo) - VNliI^o 
M = V^0a (a = 1,2, 3), 



(24) 



where we have assumed spatially homogeneous order pa- 
rameters and (•) denotes taking the operator's expecta- 
tion value. We choose 0o to be real (and positive). For a 
fixed number of bosons Nb and at T = 0, the fields obey 
J2a=Q IV'ctP = "-bj where rib = Nb/Ni, is the bosonic fill- 
ing factor (density). 

Starting from a bosonic field ^'(x) which we assume to 
be slowly varying in real-space continuum, i.e., on length 
scales larger than the lattice spacing I ai_2 1 — 1, (assuming 
low temperatures compared to Tkt) we approximate. 



^'(x) « Vo + 



(25) 



with a modulation in real-space that is solely due to 
the wavevectors Qq. We derive the free energy func- 
tional J-'b for the homogeneous system in detail in Ap- 
pendix [C] It contains the quadratic and quartic terms, 
in the {ipi}, that are invariant under all the symmetries of 
the isotropic triangular lattice. These are one 3-fold rota- 
tion, two reflection symmetries and the two translations 
by a.i [t^. Up to quartic order in the order parameters, 
it reads 



Nr 



Too 100 P + mi|i/jQp + ^ + Xlf^^i . (26) 



j=0 



(27) 



where xpg = (0i,i/'2,'03) and the different terms read 

eo = |0ol' 

ei = |tAQl' = l^i|' + l^2|^ + |^3l' 

+ 2 (iv'in^sp + i^in03p + i02n^3p) 

62 = l^oHtAgP = |0OP (IV'lP + I^2P + 
Fi = \^lf + \ij2f + \H' 

F2 = (V? + i^l + i^iy {^l + i>l + i'D 

F3 = "00 ("0102 "03 + cyclic permutations) + c.c. 

F4 = {^i + + ^ly + c.c. , 

with ■0 ~ (0Oj 0ij 02, V's)- The free energy Tb contains 
nine coefficients: two masses {mo, mi}, and seven inter- 
action parameters {uq, wi, M2, <?i, 52, 53, 54}- 

Instead of giving an exhaustive phase diagram of JFb, 
we match these coefficients with a microscopically de- 
rived bosonic mean-field Hamiltonian H^^, in the spirit 




FIG. 4: (Color online) Kagome supersolid phase on the tri- 
angular lattice. Darker lattice sites exhibit a higher bosonic 
density. The bosonic density is smaller at the lattice sites 
that belong to the Kagome-sublattice structure. If the boson- 
fermion interaction is repulsive, the fermionic density is larger 
where the bosonic density is smaller (lighter lattice sites). 



of Weiss mean- field theory [73| . We obtain H^^ from in- 
tegrating out the fermionic degrees of freedom in the full 
Hamiltonian of Eq. ([2]), as before, leading to the effec- 
tive bosonic Hamiltonian of Eq. pU|) . We then perform a 
(more general) Bogoliubov approximation of the bosonic 

operators &q = J2a=Q(^Qc)H^- Qa) +bq^Q„, where 5q 
describe the fluctuations around the mean-field values, 
and we have defined Qo = 0. Neglecting fiuctuations, 
one arrives at 



Nl 



^6(Qa)|V'o 



+ E' ^(Q^~Q^V ;0a^.0^., (28) 

Q,/3,7,(5=0 

where the second sum is restricted to Qa+Q^ = Q^+Qd". 
Since the Lindhard function possesses all the symmetries 
of the lattice, it is identical at the three nesting vectors 
q = Qi,2.3, and we can define the interaction coefficients 



(29) 



u^U{T,())^Ubb + Ulfx{T,0) 
g = U{T, Qi,2,3) = Ubb + U^fXiT, Qi,2,3) • 

This simplifies the Hamiltonian to the general form 



with 

, . , _ 61 + 292 - 2Fi + F2+ 4f3 + Fj 

^ - ^^^^ 

being a function of the direction of the vector '0/1 ■01 only. 

In this form, one easily derives the generalized crite- 
rion to avoid phase separation. Stability requires that 
the quartic coefficient is always positive, because the free 
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This demands 
0, the stabihty 



energy must be bounded from below. 
u + gW > , and because 3 > W > 
conditions are given by 

u>0, if g > 
u>3\g\, if5<0. 

Matching parameters between Eqs. and 
yields for the mass coefficients the expressions 



(a) 



mo 
mi 



6(0) 

6(Qa 



2tb 



fj-b 



(32) 

iOl) 



(33) 



We refer to Appendix |D] for the expressions of the inter- 
action coefficients. Since the chemical potential in the 
superfluid phase is given by /i^,^^' = —6tb + UbU, which 
contains a mean-field energy shift due to interactions, 
the mass = —UbU is always negative. This indicates 
that the system wants to condense into the tpo-mode in- 
dependently of the value of tt, because it does not cost 
any kinetic energy to add a boson to the zero-momentum 
condensate. In contrast, the mass mf^ = 8tb — UbU de- 
pends on the ratio of kinetic to interaction energy. It 
costs a kinetic energy amount of 8tb to add a boson into 
one of the nesting modes ■01,2,3- We will later confirm 
that for smaller hopping amplitudes tb, the supersolid 
already occurs for a smaller interaction strength Ubb- 



2. Mean-field phase diagram 

Minimization of the Hamiltonian H^^ of Eq. pO]) 
yields the phase diagram of the system. Here, we 
present only the main results, that can be obtained by 
a straightforward numerical minimization. We refer to 
Appendix [D] for a detailed analytical study. 

We numerically minimize Eq. (j30p using the most gen- 
eral ansatz ipj = r^e^'^J with rj > 0, {j — 0,1,2,3), 
4>o = 0, < 01,2.3 < 27r. We find that while for g > 0, 
the superfluid has a lower free energy than any super- 
solid phase, such a phase occurs for sufficiently negative 
g < 0, where the system tends to order in a symmetric 
way with respect to the three nesting fields. Furthermore, 
the phases of the supersolid order parameters are locked 
to the superfluid phase 0i,2,3 = 0. In this way, the 3-fold 
rotational symmetry of the system is preserved. In real 
space, the resulting density wave order has the pattern of 
the Kagome lattice, as illustrated in Fig. H) Note that this 
is also consistent with the form of the fermion-induced 
interaction between the bosons (see Appendix [B]). This 
Kagome-state can be written as 



■00 — \/nb cos 
fnb 




(34) 



with Q < 9 < tt/2. To construct the phase diagram, 
we first observe that the transition between the super- 
fluid and supersolid phase is between two ordered phases. 







SS 1 1 






SF 


PS i ^ 






9 



9 



FIG. 5: (Color online) Phase diagram from the microscop- 
ically matched bosonic mean-field theory in the parameter 
space (m, (?) for tb — O.ltf and filling factors nt — 1/4 (left) 
as well as nt = 5/4 (right). It contains regions of SS, SF 
and PS. PS occurs below the instability line u = 3\g\ (solid). 
The two vertical lines denote the critical interaction strengths 
gi^^ — —Atb/ub (dashed) and gi^^ = —2tb/nb (dotted), respec- 
tively. In between, the phase boundary occurs at gc (solid), 
which is given in Eq. ((35}. The tricritical point occurs at 

ffmax = -20tb/9nb. 



Starting from the superfluid and replacing 0o ^ in 
Eq. (l30l) . we find that the free energy contains a third 
order term in V'i,2,3- This third-order term leads to a 

(2) 

local minimum in the free energy for g < gc — —^tb/nb- 
This local minimum, however, becomes the global mini- 
mum only at the more negative value 



12tbU 



16tb + SubU 



(35) 



This equation embodies a more general condition to en- 
ter the supersolid phase than the instability criterion 
of Eq. (HH). Together with the stability conditions in 
Eq. (|32[) . it divides the parameter space {u,g) into three 
phases: superfluid (SF), supersolid (SS) and phase sepa- 
ration (PS). In Fig. [5l we show the resulting phase dia- 
gram for a fixed value of tb = O.lt/ and different bosonic 
fillings Ub = 1/4,5/4. The maximal value of g to enter 
the supersolid phase is ^max = —20tb/9nb, which is the 
tricritical point of the phase diagram. 

The superfluid-supersolid phase transition is of first- 
order [tiI, [t^] , because the curvature of the free energy 
around the superfluid minimum, i.e. the effective mass 
of the fields 0i,2,3, remains positive even below gc until 
g < gc^\ The first-order transition region between the 
two vertical lines gc'^^ — g^c^ = 2tb/nb shrinks with larger 
Ub, which is fully consistent with numerical results of Ref. 
[78| , where a first-order region could only be identified for 
sufflciently small Ub. Again, for details of the calculation, 
we refer the reader to Appendix [Dl 
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3. Transition temperatures and low-temperature phase 
diagram 

Here, we want to derive transition temperatures Tpg 
and Tgs using the mean-field criteria of Eqs. ([5^ and 
Eq. ((35)) . that generahzc the instabihty expressions of 
Eqs. (fTHl [21 ]) . This allows us to draw a low-temperature 
phase diagram of the system. 

Starting from the superfluid phase, we can calculate 
the interaction parameters u and g at temperature T 
and for fixed Ubbi Utf using Eqs. (^5)1 . Upon lowering the 
temperature, both u and g decrease and finally reach the 
superfluid phase boundary (see Fig. [5]). Phase separation 
is avoided as long as u > SigmaxI or 



x(rps,o) = -^ 



bf 



1 - 



20tb 
SribUbb 



(36) 



Note that the instability analysis only demands the less 
strict condition u > 0. The supersolid phase appears for 
5 < 5c or 



xiTssM 



1,2,3) 



-SnbUbbjUbb + cpU^f) - itbjlUbb + ScoUjf) 
Uij[mb + 3nbiUbb + coUif)] 



(37) 



where we have defined cq = x{T,0). Note that the 
instability analysis required a more negative value of 

g < g'P = -'^tb/rib. 

We extract the transition temperatures {'7ps,rss} 
from Eq. ([36|) and (I37|) . using the expressions of 
X(r, Qo,i,2,3) given in Eqs. (US]) and One should 

note that they are only valid in the superfluid phase. 
Nevertheless, they allow us to divide the (C/bb,C/&/) pa- 
rameter space into a supersolid and a phase separated 
region, in the following way. 

Assume that Tss > Tps and the system is in the su- 
perfluid phase at a temperature T > Tgg. As the tem- 
perature is lowered, it will become supersolid at T = Tgs. 
As we will show in the next Sec. IIIIBl the instability to- 
wards phase separation is removed inside the supersolid 
phase, since the fermionic spectrum acquires a gap at 
the Fermi energy. The system remains supersolid for all 
temperatures T < Tss- In contrast, if Tpg > Tss and we 
again lower the temperature starting in the superfluid 
phase at some T > Tps, the system will phase separate 
at T = Tps. Then, the fermionic density deviates locally 
from rif = 3/4 and the Fermi surface is not nested any- 
more. Hence, the instability towards supersolid forma- 
tion is removed. In Fig. [51 we show the resulting phase 
diagram and contours of constant supersolid transition 
temperatures Tss for fixed Ub — 5/4, tb = O.ltj. It is 
important to note that the transition temperatures are 
invariant under the transformation 

tb atb, Ubb aUbb, Ubf — > VctUbf . (38) 

Therefore, the same Tss can be found for other bosonic 
hopping amplitudes tb under proper rescaling of Ub / and 
Ubb- 



3.0 
2.5 
2.0 
^ 1.5 
1.0 
0.5 



SF-MI 
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I I I I I J 
1 1 1 I 



ss 



' ' ' 'y, 
^1 1 1 1 i/i 
I iiiyi / 



I". 



7 /4'/ /■ . 



-0.5 



0.0 0.5 1.0 1.5 2.0 

FIG. 6: (Color online) Low-temperature phase diagram for 
parameters {Ubf,Ubb) (in units of tf), and fixed values of 
tf — 1, tb = O.ltf, nt — 5/4, Uf = 3/4, obtained from 
tlie bosonic mean-field theory on the triangular lattice. The 
SS appears where Tss > Tps, PS appears where Tps > Tss. 
Temperatures are from Eqs. (|36p and p7p . The phase bound- 
ary (solid line) is defined by Tss — Tps- The horizontal 
dashed line denotes the critical ratio for a competing SF- 
MI phase transition, that occurs at {Ubb/tb)c ~ 26.5 [74l |. 
Other dashed lines are supersolid transition temperature con- 
tour lines corresponding (from right to left) to logjQ(rss/f/) = 
-0.5, -1,-1.5, -2,-2.25, -2.5,-2.75, -3,-4. 



B. Supersolid density wave modulation from 
fermionic mean-field theory 

So far, we have located the supersolid parameter 
regime (see Fig. [6]) and identified a highly symmet- 
ric Kagome-type supersolid phase with equal amplitude 
modulation in all three nesting wavevector modes (see 
Eq. ([23])). In this section, we use a different (fermionic) 
mean-field approach, that treats the fermions exactly, to 
calculate the effect of this condensation, i.e., '01,2,3 ^ 
0, onto the fermions. In general, one finds that the 
fermionic spectrum acquires a gap at the Fermi surface 
for nonzero ■01,2.3 • This is energetically favorable for the 
fermionic subsystem, however, notice that this energy 
gain has to be sufficient to, at least, compensate the ki- 
netic and interaction energy cost of adding a boson with 
a large wavevector in the system. 

We then determine the amplitude of the density wave 
modulation in the supersohd phase. We find that the 
amplitude is generally rather small and increases for de- 
creasing bosonic hopping amplitudes tb- For tb — O.Oli/, 
we find it to be (maximally) about Anb/rib = 0.1 at 
zero temperature and involve only about 0.1% of all 
the bosons. We like to mention that similar density 
wave modulations were predicted in Refs. jssl . pTSj using 
Dynamical Mean-Field theory (DMFT). As we show in 
Sec. IV C[ the experimental detection of such a small den- 
sity wave is not feasible with current technology, however, 
we note already at this point, that one finds significantly 
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larger density wave amplitudes for the square lattice ge- 
ometry (see Sec. lIVp. 

We begin by replacing the bosonic operators 6qo, 1.2,3 
with ^/'o,i,2,3 as in Eq. (|24p. Neglecting any fluctu- 
ations, the full Hamiltonian of Eq. ([2]) then becomes 
Hf = nf^ + nf^ + Hf \ where 



0.20 



H 



(1) 3 



a, /3, 7, (5=0 



(39) 



qG-BZ 

where the primed sum is restricted to Qc -I- = + 
Qi, and we have incorporated a mean- field energy shift 
of the fermions due to the presence of the bosons given 
by C/b/«b X)q /q/q iiito the chemical potential /i/. 

The first term H^'p describes the kinetic and inter- 
action energy of the condensed bosons. The kinetic en- 
ergy cost to take a boson from the superfluid condensate 
ipQ and add it to one of the nesting modes is given by 
C&(Qi,2,3) = 8<b. The (quartic) interaction energy ob- 
viously increases with the number of nonzero fields "0^ 
and larger Ubb- To find out whether these energy costs 
of having nonzero '01,2,3 can be compensated by the last 
two (fermionic) terms, we diagonalize them. 

For this, we need to symmetrize the expressions 

with respect to adding a nesting vector Qq, such 

(2) (3) 

that, in matrix form, we write: H)' + Hf' = 
ELEa,/3/k+Q„^a/3/k+Q^, where h^p is given by 



hafS = Sai3^f{k + Qa) 

+ (1 - S^p)Ubf (C# + r^^s + c.c.) , (40) 

where the sum over wavevectors is restricted to 1/4 of the 
first Brillouin zone, and (7, S) ^ (a, (3) G {0, .., 3} as well 
as 7 7^ (5. By diagonalizing /iq,^, we obtain the fermionic 
eigenenergies S(k, {ipa})- 

The free energy at finite temperatures T > reads 



^ = E^''(Q")iv^"i' + ^ E' r^r.i^yi^s 



T 

iv7 



Ein(i 



E(k,{^„})/T 



) 



(41) 



As the temperature goes to zero (T 0), the last 
term becomes a sum over the lowest Ni^Uf eigenstates: 

l^EiIwr/tj=i^(^'{V'a}). Since the supersolid phase 
appears at temperatures that are much smaller than the 
fermionic hopping amplitude, Tss tf, results of finite 
and zero temperature calculations agree. 

We calculate T at zero temperature on a finite lattice 
with Nl = 500 x 500 sites, as a function of the {0q,}, 




O.OQ 



0.1 0.2 0.3 0.4 0.5 0.6 
Ubf 

FIG. 7: (Color online) Density wave modulation Ant/nt as 
a function of Ubf (in units of tf) for fixed Utt ~ 0.25t/, 
rii, = 5/4, ti, = O.Oltf. The vertical dashed line denotes the 
phase boundary to the phase separated regime for this value 
of Ubb- 



using the Kagome ansatz that we have used before: 0o = 
y^cosB, 01,2,3 = y^e*"^ sin6', with angles < 9 < 
7r/2 and < < 27r. We then determine its minimum 
as a function of {9, 0}, for a certain choice of parameters 
{Ubb, Ubf,tb,nb}. 

We find that the minimum always occurs for = 0. 
The location of the minimum as a function of 9 deter- 
mines whether the system is superfluid or supersolid. If 
the minimum occurs at 9 = 0, condensation into the nest- 
ing modes is energetically not favorable and the system 
remains superfluid. In contrast, if it occurs at 9 > the 
system can lower its energy by establishing (6q^) ^ 
and becomes supersolid. The value of 9 determines the 
amplitude of the density wave in the supersolid phase. 

At zero temperature, the minimum occurs at 9 > for 
any parameter pair {Ubb, Ubf) in the weak-coupling region 
(see Fig.|6]). However, the amplitude of the corresponding 
density wave modulation in the supersolid phase, which 
we define as 



Aub = maxi{blbi) - mmi{blbi 



(42) 



with the bosonic operators approximately given by bi ~ 
-00 -t- X]a=i ^aS.^^" '^^ , varies significantly with the values 
of Ubb and Ubf- 

More precisely, Artf, is proportional to Tgs/Ubf, a 
relation that can be understood from Bardeen-Cooper- 
Schrieffer theory, where the gap at zero temperature is 
proportional to the superconducting transition tempera- 
ture. Here, the gap in the fermionic spectrum is deter- 
mined by the product Ubfip'^tpis with a ^ 13 G {0, ..,3}, 
as can be seen from Eq. (|40)) . For the Kagome supersolid 
we find, 



Anb 8 

Ub 



■ sm I 



(V3 



cos ( 



sm I 



(43) 



which takes values between < Aub/ub < 4 and is shown 
in Fig. [71 For parameters Ub — 5/4 and tb — O.Oli/, 
it takes values up to Aub/ub = 0.1. We also ex- 
tract the ratio AubUbf /Tgs ~ 9 ± 1 using Tss from 
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(a) 




FIG. 8: (Color online) Striped supersolid phases on the square 
(a) and anisotropic triangular (b) lattice. The density wave 
is characterized by a single nesting vector: (a) Qsq — (vr,7r), 
and (b) Q3 = [0,2^/^/3). Like in Fig. |4l darker lattice sites 
exhibit a larger bosonic density. 



Eq. p7p . For comparison, we mention that one can de- 
rive an exact analytical relation on the square lattice [2l[ : 
AubUbf/Tss = 4/Ci « 3.53 (see Eq. where Tgs is 

defined in Eq. dST 



2d anisotropic square lattice, however, the presence of 
divergencies close by in energy, lead to a significantly in- 
creased value of the density of states at the Fermi surface 
energy (see Fig. [^ . 

In general, compared to the triangular lattice, the 
square lattice shows supersolidity at higher temperatures 
Tgs, and with larger density wave modulations Arib, be- 
cause the nesting relation is fulfilled for more k-values. 
Therefore, the square lattice geometry is experimentally 
advantageous over the triangular one, and will be dis- 
cussed first. 



A. Square Lattice 

An instability analysis and a fermionic mean-field the- 
ory of the supersolid phase on the isotropic square lat- 
tice has been discussed in Ref. 



21| . Here, we will in- 



clude a bosonic Landau-Ginzburg mean-field treatment 
and also generalize to the case of spatially anisotropic 
hopping amplitudes. To obtain a clear presenetation of 
the main results, the calculational details can be found 
in Appendix [El 



IV. STRIPED SUPERSOLID PHASES ON 
ANISOTROPIC LATTICES 



1. Instability analysis 



In this section, we discuss the case of anisotropic hop- 
ping amplitudes on the two-dimensional triangular and 
square lattices, where only one nesting vector remains as 
a result of the reduced symmetry. The isotropic square 
lattice also exhibits a single nesting vector. As a result, 
we find supersolid phases that show a striped pattern in 
real-space; see Fig. [5] 

On the triangular lattice, there are two ways to in- 
troduce anisotropic hopping (see also Fig.l^a)). Either, 
the hopping amplitude tfi > tf2, leading in the limit 
tfi 3> tf2 to an array of weakly coupled chains, or the 
opposite case of tfi < tf2, where the system resembles 
the isotropic square lattice in the limit tfi ^ tf2. We 
will therefore discuss the square lattice geometry in some 
detail in Sec. HV^ 

Whereas for the triangular lattice, nesting is always ac- 
companied by the occurrence of van Hove singularities (at 
the same fermionic filling Uf), and thus x(T, q) diverges 
at zero and the nesting vector, one can separate both 
phenomena to occur at different n/ on the anisotropic 
square lattice (compare Figs. [51 and [T^. As a result, on 
the anisotropic square lattice, the Lindhard function is 
regular at q = even in the presence of a nested Fermi 
surface, and the tendency towards phase separation at 
low-temperatures is suppressed. On the other hand, we 
will find that the divergence at the nesting vector be- 
comes weaker, which leads to slightly smaller supersolid 
transition temperatures Tgs- We note that this situation 
is similar to the case of three dimensional lattices. There, 
critical points on the Fermi surface are integrable and the 
density of states is regular everywhere. In the case of the 



The fermionic (bosonic) dispersion relation for the 
square lattice is given by 

C/(b)(q) = -l^f{b) - 2[t/(f,)i cosqi -I- t/(b)2 cos 92] ■ (44) 

It is useful to introduce the anisotropy parameters ryj-^-) = 
tf(b)i/tf(jb)2 and work with a dimensionless chemical po- 
tential /2/ = iif/2tf2. 

The Fermi surface for different values of and jlf 
is shown in Fig. [51 There are two special values of the 
chemical potential: first, jlf — Q (dashed contour), where 
the system is particle-hole symmetric and shows perfect 
nesting at wavevector Q^g = (tt, tt) for any value of r/. 
Second, iov jlf = ±(1 — r/), the density of states has a van 
Hove singularity due to the critial points q = (±7r, 0) on 
the Fermi surface (dotted contour). Only in the isotropic 
system, these values are equal and given by /x/ = 0. 

The density of states g{jLf,rf) for the isotropic and 
anisotropic lattice exhibits a logarithmic singularity at 
jlf — ±{1 — rf). It can be calculated analytically (see 
Appendix IE 1 ap . which also contains a plot of g{jif,rf)). 
The anisotropic density of states is regular at jif — 0, 
where it is equal to 



g{jlf = 0,rf<l) = 2NoKirf)., 



(45) 



where Nq = l/27r^t/2 and K{x) denotes the complete 
elliptic integral of the first kind. 

Next, we calculate the Lindhard function x{T, q) 
(Eq. On the isotropic square lattice, it was shown 

in Ref. 21] that the Lindhard function diverges at q = 



12 



0; Qsq &S 



x(T,o) 



2 



16Cii 



/ 



T 

16Cii 



In 



n 2 



(46) 
(47) 



due to the combination of van Hove singularities and 
nesting. We note that as opposed to the triangular lat- 
tice, here, the minimum of the Lindhard function close 
to the nesting vector always occurs at q = Qsg (see 
Sec.|llE|. 

On the anisotropic lattice, the density of states is reg- 
ular at p,f — 0, and therefore the Lindhard function 
is regular at q = 0: x{T,0) = —g{0,rf). The ten- 
dency towards phase separation at low temperatures is 
removed for rf < 1. On the other hand, the divergence 
of x(r, Qsq) in the absence of a van Hove singularity at 
p,f = 0, is only linearly logarithmic: 



x(T,o) = 



^^■'^(^^^ 



-2NoK{rf) (48) 



deg{e) 



tanh(e/2r) 



-2iVo/-i:(r/)ln 



-2e 

4e^{l + rf)h{rf)tf2 



ttT 



(49) 



The fitting function h{rf) = oq + airf with ao = 1.96, 
fli = —1.67 occurs from comparing numerical results to 
an analytical approximation that replaces g{e) ~ g(0) 
in Eq. (|49p . which neglects the (nearby) divergence in 
the density of states at /2/ — ±{1 — rj). The slope of 
x(T, Qsq) is given by the regular density of states at the 
Fermi surface. 

These divergences result in two instabilities, as we have 
shown in Sec. [iTl For the isotropic lattice, phase separa- 
tion occurs at the temperature, 



^^s"(^/ = 1) = 16Cit/exp 



Ubb 



whereas the supersolid transition temperature reads 



Tssirf = 1) = 16Cii/exp 



/ 2Ubt 



44 \ 
nbUbb) 



(51) 

with Ci = 2expC/7r ~ 1.13. 

For the anisotropic lattice, phase separation only oc- 
curs above a critical interspecie interaction strength, 
which is given by 



PS,inst. 



Ubb 



(52) 



The supersolid transition on the anisotropic lattice occurs 




FIG. 9: (Color online) Fermi surfaces of the isotropic (a) and 
anisotropic (b) square lattice with anisotropy parameter r/ = 
0.75, for different fillings. Nesting occurs at /i/ = (dashed 
lines) and van Hove singularities at /i/ = ±(1 — r/) (dotted 
lines) . 



at a temperature 
Tss{rf < 1) 
= At f2 exp 



-Ubb 



2(^61 + tb2) 

mUbb 



(53) 



9iO,rf)U^f 
where A = 2Ci(l + rf)h{rf). 

2. Phase diagram from bosonic mean-field analysis 

The effective bosonic mean-field Hamiltonian on the 
square lattice (see Sec. IIII Al) reads 



rroff 

Nl 



1 

E 

a=0 



1 



UQc.)\^c.r + -{u + gW.,)\'iP\ 



(54) 



with u = U{T, 0), 5 = U{T, Qsq) and 



(V-oV^l* + C.C.)2 



sin2(26')cos2(0), (55) 



where we have parametrized the bosonic fields as xp = 
(50) (V'o,V'i) = (V^cos6',y?i^e*'^sin6') and Qo = 0, Qi 

Qsq. For arbitrary (j), we can distinguish between three 
different phases: the cases 9 = 0, 7r/2 correspond to a 
superfluid and a pure density wave phase, respectively. 
For < 9 < 7r/2, the system is supersolid. Furthermore, 
the stability requirement, which is that the fourth order 
term should be bounded from below, is given by u > 
for g > 0, and u> —g ior g < 0, because < Wsq < 1. 

To obtain the phase diagram, we minimize the Hamil- 
tonian. This is done analytically in Appendix IE 1 bl 
where we find that the pure density wave phase always 
has larger energy than the superfluid, which is the ground 
state for g > 0. Comparing the energy of the superfluid 
with the supersolid, we find that the supersolid occurs for 
sufficiently attractive interactions at 



9 < 9sq,. 



2(<61 + tb2) 



rib 



(56) 
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FIG. 10: (Color online) Phase diagram for the isotropic 
square lattice as a function of {Ubb,Ubf) (in units of tf) 
with fixed tt — 0.39t/, ni, = 3/2, Uf = 1/2, so that 
fif — 0. Solid and dotted lines indicate the phase bound- 
ary obtained from the bosonic mean-field calculation and 
instability analysis, respectively. They separate a SS from 
PS regions. Horizontal dashed line denotes the critical SF- 
MI ratio {Ubb/tb)c = 16.5 [t^. Other dashed lines indicate 
constant supersolid transition temperatures: logiQ(Tss/t/) — 
0, -0.25, -0.5, -1, -2 (right to left). 



In contrast to the triangular lattice, the superfluid- 
supersolid phase transition on the square lattice is of 
second order. The resulting supersolid transition temper- 
ature thus coincides with the transition temperature ob- 
tained by the instability analysis Tss of Eqs. ([5T|) and ([55]) 
for both isotropic and anisotropic lattices. 

However, the condition to avoid phase separation is 
modified from it > to u > —gsq.c, leading to a larger 
transition temperature towards phase separation on the 
isotropic lattice 



Tpg = 16Cit/exp 



Ubb 



I - 



44 



(57) 



which obeys Tps > T^^'i^. 

On the anisotropic lattice, the critical interspecie in- 
teraction strength to avoid phase separation is modified 
as well, and reads 



rrPS _ 



' Ubb - 2(4i +tfc2)/wfc 
2NoK{rf) 



(58) 



The resulting phase diagrams for the isotropic and 
anisotropic square lattice as a function of {Ubb, Ubf}, to- 
gether with contour lines of constant Tgg, are shown in 
Figs. [To] and [ni 



Density wave modulation from fermionic mean-field 
theory 



Following the analysis of Sec. IIIIBl we calculate the 
amplitude of the supersolid density wave by a fermionic 
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FIG. 11: (Color online) Phase diagram for the anisotropic 
square lattice with r/ = ri, = 0.75, as a function of {Ubb, Ubf) 
(in units of tf) and with fixed tb2 ~ 0.39f/2, nb — 3/2, 
fif — 0. Solid and dotted lines indicate the phase bound- 
ary between SS and PS obtained from the bosonic mean-field 
calculation and instability analysis, respectively. The hori- 
zontal dashed line denotes the critical ratio {Ubb/tbijc = 16.5 
to the competing MI phase. Other dashed lines indicate 
constant supersolid transition temperatures log]^Q(rss/t/) = 
-0.1,-0.25,-0.5,-0.75,-1,-1.5,-2.0 (dashed, from right 
to left). 



mean-field approach. We replace the bosonic operators 
^o.Qi by the complex fields -00 — and -01 — exp(i0) 
(ro,i € M) and diagonalize the resulting fermionic spec- 
trum exactly. This way, we derive the finite temperature 
free energy 



TLbUlf 



-T 



E 

k,s = ± 



In 1 



C/febA^ cos^ ( 

-H(k,A),/T 



(59) 



which contains the gap A = 2Ubfrori and the fermionic 
eigenenergies 



S(k,A)± =±./^/(k)2-|-A2cos2, 



(60) 



Details of the derivation are given in Appendix IE 1 cl 
From d^J- ~ 0, one finds cf) — rmr with integer m, and 
the gap equation arises from d^J^ = as 



2 + TB 

A 



BF 



2 



k 



tanh[S(k, A)+/2T] 
S(k:A)7 



(61) 



where Ts = 4(^{,l+^^,2)/?^6^76b, >^bf = NoU^j/Ubb and the 
summation is over 1/2 of the first Brillouin zone. Solv- 
ing for the supersolid transition temperature by setting 
A (Tss) — 0, reproduces the results from the instability 
analysis and the bosonic mean-field theory (Eq. ([5T1 [55)1 . 

The density modulation in the supersolid phase Ant 
is proportional to the gap 



2A 

Arib = — , 
Ubf 



(62) 
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FIG. 12: (Color online) Normalized density wave amplitude 
Anb/rib = (max{6j&i) — mm{blbi))/nb in the supersolid as 
a function of Ubf (in units of tf2) for (a): isotropic square 
lattice with tb = 0.39 1/, Ubb = 5.7 tf, Ub = 3/2, and (b): 
anisotropic square lattice with r/ = 0.75, tb\ — 0.06 1/2, tb2 ~ 
0.09f/2, Ubb ~ 0.9t/2, rib = 3/2. This choice of parameters 
corresponds to a particular experimental realization that will 
be discussed in Sec.|VB](see Tab. Iml |IV|| . 




FIG. 13: (Color online) Fermi surfaces of the anisotropic tri- 
angular lattice for the hopping amplitudes tfi = 1, tf2 ~ 0.75 
(a) and t/i = 1, t/2 = 1.25 (b) and different chemical po- 
tentials fif. Nesting occurs for fif — 2tfi (dashed) with a 
single nesting wavevector Q3 = (0, 27r/\/3)- It is accompa- 
nied by van Hove singularities in the density of states due to 
critical points at q = ±Qi,2- Critical points also occur for 
fj,f /tf2 = 4: — 2rf (dotted). The hexagon (thick line) denotes 
the first Brillouin zone. 



because the expectation value of the number operator 
in the supersolid state reads = |'0oP + IV'iP + 

cos(Qi • Xi). We can solve for the zero temperature 
gap A(T = 0), using the modified density of states in the 
gapped system G(S, A) = ^(VS^ - A2)|S|/\/S2 - A^. 
For the isotropic lattice, one finds A(0) = 2Tss{rf = 
1)/Ci (see Eq. ([5T|)). and in the case of the anisotropic 
lattice, we find A(0,r/) = V2Tss(r/ < l)/Ci/i(r/) (see 
Eq. ([55)1 ). The density modulations at T = thus read 
for the isotropic lattice 



Aub 



4Tss(r/ = 1) 



and for the anisotropic lattice 

2V2Tss(r/ < 1) 



Aub = 



CMrfWbf 



(63) 



(64) 



where, as defined earlier, h(rf) — 1.96 — 1.67r/. In 
Fig. I12[ we show typical density wave amplitudes Anb 
as a function of Ubf for fixed Ubb and tb- 



of van Hove singularities at the same chemical poten- 
tial fj,f. Therefore, also x(T, 0) diverges logarithmically 
at low temperatures and supersolid formation competes 
with phase separation. 

To show that only -00,3 tend to condense, we derive 
the effective bosonic mean-field Hamiltonian 



(65) 



+ -{u + vV + gW)\ip\'^ , 



with the masses rrii = ^^(Qi) and interaction coefficients 
u = t/(r, 0), V = UiT, Qi,2), 9 = UiT, Qg), and 

+ [■0oV'lV'2V'3 + ^oV'l'02V'3 + 2V'oV'tV'2V'3 + C.C.] 

^ {Ms + i^ir2 + c.c.f . 

(66) 



B. Triangular lattice 

For spatially anisotropic hopping on the triangular lat- 
tice, the Fermi surface is still nested at the chemical po- 
tential fif — 2tfi. As shown in Fig. 1131 however, only 
one nesting vector Q3 = (0,27r/'\/3) remains. This is 
true for both i/i < tf2 and tfi > tf2. As a result, the 
Lindhard function only diverges for wavevectors close to 
q = Q3 (and no longer at q = Qi,2)- We show below 
that the tendency to condense into the modes ipi^2 is then 
removed. With only i/jq.s being nonzero, the supersolid 
has a striped pattern in real space; see Fig. [51 In addi- 
tion, nesting is always accompanied with the occurrence 



We analyze the Hamiltonian in detail in Appcndix lE 21 
where we find that, again, kinetic energy considerations 
will select the superfluid state for positive g. For negative 
g, the system can possibly lower its energy, compared 
to the superfluid, by allowing for nonzero ■03 while still 
having = -02 = 0. In this subspace of possible field 
values, the mean- field Hamiltonian is of the same form 
as on the square lattice and we refer to the discussion in 
Sec. IIVA2I and Appendix [ETbl 

Numerical evaluation of the Lindhard function shows 
that, for r/ < 1 (r/ > 1), the supersoHd transition tem- 
peratures are above (below) the ones for the isotropic 
triangular lattice. They are always smaller than Tgs on 
the isotropic square lattice. 
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C. Summary 

On the square lattice, we have found a checkerboard- 
type supersolid that occurs on the isotropic lattice at 
temperatures as large as Tss — tf = Tf/4: (for ni, = 3/2). 
For the anisotropic square lattice, we have found only 
slightly smaller temperatures Tss — 0.6 ~ Tp/h. 
There, phase separation only occurs above a critical in- 
teraction strength U^^. 

The fermionic energy gap A, which is proportional to 
Tgs, is related to the supersolid density wave via Arif, = 
2A/[/t,/- Hence, larger amplitudes occur for smaller 
bosonic hopping, since, according to the transformation 
of Eq. ([55]) , the same Tss then occurs at a smaller Ub / . 
The amplitudes are similar for isotropic and anisotropic 
lattices. For ti, = 0.39^/, we find Arib/ni, « 0.5, and for 
the smaller value of ti, « O.lt/, we find Aub/ni, « 0.9. 

On the anisotropic triangular lattice, we have iden- 
tified a striped supersolid that competes with phase 
separation. For < l(r/ > 1), it occurs at larger 
(smaller) temperatures than the Kagome-supersolid on 
the isotropic lattice. 

In order to compare different lattice geometries, we es- 
timate an upper bound for the supersolid transition tem- 
peratures as the temperature Tgg, where 

x(rs*s,o) = x(rs*s,Q.)., (67) 

where i = 1(3) for the square (triangular) lattice. 

By numerically computing xC^j i) o^i the triangular 
lattice for different anisotropy parameters rf — i/i/t/2, 
we observe that Tgg increases with decreasing values of 
rf < 1. For the isotropic triangular lattice, one finds 
Tgg(r/ = 1) « 0.2 1/2- For the isotropic square lat- 
tice, which is the limiting case of ^ 0, one calculates 
Tgg(r/ 0) K 1.2tf. In between, one has 

0.2<Zis(l^II^<1.2. (68) 

In contrast, for r/ > 1, the temperature Tgg decreases. 
I.e. T^sirf > 1) < Tisivf = 1). 

The supersolid transition temperature increases with 
n-b, and with Ubf being close to the phase bound- 
ary of supersoHd to phase separation and Ubb/h close 
to the superfluid to Mott-transition ratio C/hb/tblsF-Mi- 
However, the weak-coupling requirement \bf/tb — 
UbMoU^f/Stb < 1 obviously restricts the maximum value 
of nfc. We find that the transition temperatures consis- 
tent with Xbf/tb < 1 are close to the upper bound Tgg 
for the square lattice, but generally much smaller for the 
triangular lattice. The difference there, is typically an 
order of magnitude. 

In conclusion, the isotropic square lattice exhibits the 
highest supersolid transition temperatures Tgg of all the 
geometries considered here, with Tgg — tf = Tf/4. 



Ao[nm] 7o/27r[MHz] J,at[mW/cm"] 





670.96 


5.92 


2.56 




589.16 


10.01 


6.40 


«K 


766.70 


6.09 


1.77 




780.24 


5.98 


1.64 



TABLE I; Atomic properties of ^^Na, "^Li, **''Rb, ""K. Tran- 
sition wavelength Ao, natural linewidth 70 and saturation in- 
tensity Jaat determines ratio Vt/Vf. [S^ 

V. EXPERIMENTAL PREDICTIONS FOR 

»^Rb«K AND ^-'Na'^Li 

In this section, we will present results for specific 
experimental realizations of Bose-Fermi mixtures on 
the isotropic triangular as well as the isotropic and 
anisotropic square lattices. We predict the supersolid 
parameter regime for a mixture of ^'^Na^Li [s^, [8l| and 
of s^Rb^'OR [sl, m Hi- We also show that the unam- 
biguous experimental detection of the supersolid phase 
via time-of-flight imaging (TOF) is feasible for the square 
lattice geometry. Additional coherence peaks at the nest- 
ing vector occur with a size of up to « 0.02 (measured 
relative to the main superfluid peak) for both mixtures. 
Since the weight of these coherent peaks is reduced to 
about < 5 X 10~^ in the triangular lattice case, it proves 
more challenging to detect supersolidity in this case. We 
therefore propose a combination of usual time-of-flight 
absorption imaging with noise correlation techniques [26| 
to reveal the supersolid phase. 

A. Relating Hamiltonian to experimental 
quantities 

The Hamiltonian of Eq. ([1]) contains the parameters 

C/fcb, C/fc/,/i/,/ib} , (69) 

which can be expressed by the microscopic and experi- 
mentally tunable parameters 

{to/ , TOfc , rt/ , Tifc , abb ,abf,X, Vf , ,Vf} . (70) 

Here, TOj(f,) is the fermionic (bosonic) mass, n/(6) is the 
fermionic (bosonic) density, and abb(bf) the s-wave scat- 
tering length of boson-boson (boson- fermion) interaction, 
that can be tuned by an externally app lied static mag- 
netic field via Feshbach resonances [86[. A is the wave- 
length of the optical lattice laser, Vjj^'^ the optical lattice 
laser intensity in the x, y, z-direction, given in units of the 
fermionic/bosonic recoil energy -EJ/^ = / {'^'mf /b^?), 
respectively. Here, we focus on the rectangular geometry 
for notational convenience, but the generalization to the 
triangular geometry is straightforward. 

The two-dimensional setup is realized by strongly 
quenching inter-plane hopping via VLj^-, ^ ^f(b)^ ^"^"^ ^"-"^ 
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Species A[nm] Vt/Vf 



^i, dbb [ao\ ai,/[aoj 





1064 


1.9 


1.41 


368 nK 


62 


13 




755 


fti 1 


420 iiK 


193 nK 


100 


-284 




1064 


« 2.5 


211 nK 


97 nK 







TABLE II: Ratio of optical lattice potential for bosons 
and fermions Vb[E^]/Vf[E^], measured in units of the re- 



spective 



recoil energies Elj, and coUisional properties of 
Na^Li |8l|,[8l|ll[^, ®''Rb*°K O,!!, 86]. {att,abf) denote 
the scattering lengths away from any Feshbach resonance (in 
units of the Bohr radius ao). 



f/b- 



The 



the isotropic lattice one sets V^^f^ 
lattice constant is given by A/2. 

The ratio of lattice depths experienced by bosons and 
fermions, respectively, is determined by 



Vf ^ [A 



[A-i-Ao(/)-i]7o(&)/sat(/) £^/ 
-Ao(6)-i] 7o(/) /sat(&) El 



(71) 



where Xo{f /b) is the wavelength of the relevant 
fermionic/bosonic transition, 7o(//^) its natural 
linewidth and Isatif/b) its saturation intensity. The 
prefactor C is of order unity and determined by a ratio of 
Clebsch- Gordon coefficients of the relevant transitions. 
The experimental values for the relevant transitions in 
^SNa'^Li and ^'^Rh'^^K are given in Table HI 

The resulting ratios Vb/Vf, together with the recoil en- 
ergies and collisional properties can be found in Table HTl 
Note that one can tune either one of ai,{,, ai,f over a wide 
range by applying an external magnetic field close to a 
Feshbach resonance (86| . In the following, we will con- 
sider the case where abf is tuned, leaving abb fixed to the 
(off- resonance) value given in Table |TT1 

The fermionic (bosonic) hopping amplitude in a direc- 
tion where the optical lattice depth is given by V^/(f)), can 
be expressed in closed form, if the Wannier state of the 
lowest Bloch band is approximated by a Gaussian 



However, since this approximation fails for Vj 



(72) 



fib) 



< 



lOE'^^i^y we calculate the hopping amplitudes from the 
width of the lowest energy band W{Vf(^b)) via 

tm^w{Vf(^b))/^- 

The two interaction parameters are given by [l] 
Ubb 



(73) 



Ubl 



4V2^ ^(T/-^^?^T4-)i/4 
A 



(1 + /fw^)'/' ^ ^ 



(74) 



The fermionic (bosonic) chemical potential At/(b) is de- 
termined by the number of fermions (bosons) in the sys- 
tem N^nf^b)- If the system is exposed to an overall har- 
monic confinement as is often the case experimentally. 



T (units of Tf) 
0.008 




30 32 34 36 38 40 
ay (units of ao) 

FIG. 14: (Color online) Mixture of ''°K**''Rb on the triangu- 
lar lattice: transition temperatures (in units of Tf = 6tf) 
towards supersolid Tss (solid) and phase separation Tps 
(dashed) as a function of scattering length atf (in units 
of the Bohr radius ao) Vertical dashed line indicates phase 
boundary between the supersolid and phase separation. For 
T > T'ps,Tss, system is superfluid (SF). Bosonic filling is 
Ub = 3.25 and other parameters are A — 755 nm, Vf = 7.5, 



f 



20. Inset shows abf ^ —^bf symmetry. 



the chemical potential depends on the spatial location, 
which can be dealt within the local density approxima- 
tion (LDA). For simplicity, we restrict ourselves to the 
homogeneous case, which can, in principle, be realized 
experimentally by compensating the overall confinement 
by a blue detuned optical lattice [s^l or by working with 
an external box-like potential. 



B. Experimental phase diagrams 

If we choose a particular mixture, the wavelength of 
the optical lattice A and an external magnetic field value 
that is far away from any Feshbach resonance of abb, there 
remain only (Vp^ , nb,abf) as free parameters. From the 
phase diagrams in Figs.[6l[Tni[IIl we know that the max- 
imal values of Tss are to be found where Ubf is close 
to the supersolid-phase separation phase boundary and 
the ratio Ubb/tb ~ Ubb/tb\sF-ui- Therefore, we determine 
Vj''^ (for a certain choice of Ub) by maximizing the ratio 

Ubb/tb under the constraints Ubb/tb < iUbb/tb)\sF-Mi and 
Xbf/tb < 1 (weak-coupling). 

The finite temperature phase diagram for a 
mixture as a function of the remaining free parameter 
abf is shown in Fig. [14] for the triangular lattice and in 
Fig- ESI for the isotropic and anisotropic square lattices. 

We have normalized the temperature scale by the 
Fermi temperature of the lattice Tp ~ tf, and find a 
maximal value of Tgs/Tp — 0.004 for the triangular and 
Tss/Tf — 0.2 for the square lattices. The large differ- 
ence in transition temperatures reflects the fact that the 
low-temperature divergence of the Lindhard function at 
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FIG. 15: (Color online) Mixture of '"'K^'^Rb on the isotropic 
(a) and anisotropic (b) square lattice. Shown are the tran- 
sition temperatures towards supersolid Tss (solid) and phase 
separation Tps (dashed, only in (a)) as a function of scat- 
tering length at/. Phase boundary between supersolid and 
phase separation is denoted by the vertical dashed line, (a); 
Parameters are n\, — 3/2, A — 755nm,V/ = 7.5, Vf = 20. 
Temperature is in units of Tp = 4t/. Inset shows atf ^CLtf 
symmetry, (b): Anisotropy parameter r/ = tf-i/tf2 ~ 0.75 
realized by lattice strengths Vfi — 7.5, V/2 = 6.4. Other pa- 
rameters are the same as in (a). Temperature is in units 
of Tf ~ 2(i/i + tf2)- Phase boundary occurs at (7^7 = 
y'lUbi - 2{hj+h2)/nt]/[2NoK{rf)] [see Eq. JSIJ]. 



the nesting vectors x{T, Qi) is weaker for the triangular 
than for the square lattice. We note that an estimate of 
an upper bound of Tgs is given in Sec. IIV CI 

In Tabic [1111 and IIVI we summarize the optimal choice 
of experimental parameters V/, af,/, which corresponds to 
the highest supersolid transition temperature Tss for dif- 
ferent mixtures, optical lattice wavelenghts A and bosonic 
fillings rif), for the cases of isotropic triangular and square 
(Tab. HIH) and anisotropic square lattice (Tab.[lV|. 



C. Detection of supersolid phase 

The supersolid phase can be detected unambiguously 
by time-of- flight absorption imaging (TOF), where atoms 
are suddenly released from the trap and expand ap- 
proximately freely. After a certain expansion time s, 
the observed spatial density distribution of bosons, av- 
eraged over several images, (A/'j (x)), is proportional to 
the momentum distribution in the lattice [l3| : (A/'^" (x)) = 
(TOh//is)^|?i'(k)pp(k), where 



27r 



(75) 



is the Fourier transform of the Wannier function of the 
lowest Bloch band, and the (dimensionless) momentum 
k is related to the spatial position in the cloud by k = 
Amfcx/2?is. The Fourier transform of the one-particle 



density matrix 
P(k) 



1 



'i(xj-Xfc)-k 



{b%) 



(76) 



measures the first-order coherence properties of the sys- 
tem. Here, Xj is the (dimensionless) vector to lattice site 
j, i.e. = jiai -I- j2a2 with |ai| = 1. 

For the supersolid phase, one finds that the bosonic op- 



erators can be approximated by bj ~ t/jq + J2a '^'^ 



(see Eq. ()25p . and the normalized momentum distribu- 
tion p(k) = p(k)|w(k)|VA^L|w(0)|2 takes the form 



P(k) 



^ g-|k|V[^"^A7: 



k,G„ 



(77) 

We see that a nonzero value of the bosonic density wave 
field V'q gives rise to additional coherence peaks at the 
nesting vector Qq,, where a — 1,2,3 for the isotropic tri- 
angular lattice and a = 1 for the anisotropic triangular 
and the square lattice. On the other hand, the superfluid 
component "00 manifests itself by peaks at the recipro- 
cal lattice vectors Gm = m-iGi + m2G2, with integer 
m = (mi, 777.2). Note, that this includes = 0, and 
the reciprocal basis vectors read Gi = (27r, 0), G2 = 
(0, 27r) for the square lattice, and Gi — 27r(l, l/\/3), 
G2 = 27r(— 1, I/a/S) for the triangular lattice. 

The number of atoms in the peaks at is propor- 
tional to iV'aP, which can be expressed by the amplitude 
of the density modulations An;, (see Eqs. (l43l l63l [64)) ). 
For the Kagome supersolid on the triangular lattice (see 
Eq. p4|) ). one can approximate 



1,2,3| 



T 



4n6 



(78) 



which is valid if one can neglect quadratic terms in 0, 
i.e. IV'oP ~ nb in Eq. ([34| . For the square lattice, where 
< Arib < 2n6, one finds 



^ V ^ ( 2nb ) 



(79) 



Whether the supersolid peak can be detected experi- 
mentally, is determined by its weight relative to the su- 
perfluid peak 



_ p(Q„) _ IV-apexphi^] 



IV'O 



(80) 



We therefore include this ratio in Tables IIIII and IIVI For 
comparison, we give the size of the first higher order 
superfluid peaks at the reciprocal lattice vectors. For 
a lattice depth of 14 = 11.4, one finds for the square 
lattice p(Gi,2)/p(0) = 0.31 and p(Gi + G2)/p(0) = 
0.09. For the triangular lattice, their size is given by 
p(Gi,2)/p(0) = 0.21, and p(2Gi,2)M0) = 2 • 10-^. 
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Triangular Lattice 



Species 


A[nm] 






fit 




Anb/rib 


Pa /pO 


h/tf 


Ubt/tf 


Ubb/tb 


Ubf/tf 


^Bp/rB 


«^Rb«K 


755 


9.0 
7.5 


26 
38 


1.25 
3.25 


1 X 10~^ 
4 X 10"^ 


0.01 
0.02 


1 X 10"'' 

4 X 10"* 


0.38 
0.39 


9.0 
5.7 


23.6 
14.6 


3.6 
3.3 


0.40 
0.88 




1064 


4.2 

3.5 


72 

95 


1.25 
3.25 


1 X 10"^ 
4 X 10"^ 


0.03 
0.04 


9 X 10"'' 
1.6 X 10"^ 


0.10 
0.12 


2.4 
1.8 


23.8 
14.5 


1.9 
1.9 


0.42 
0.87 


23Na«Li 


1064 


6.5 

5.8 


32 
42 


1.25 
3.25 


1 X 10"^ 
5 X 10"^ 


0.03 
0.07 


1.0 X 10"^ 

5.1 X 10"^ 


0.07 
0.08 


1.6 
1.2 


23.1 
16.3 


1.5 
1.6 


0.40 
0.99 


Isotropic Square Lattice 


Species 


A[nm] 


Vf 


abf 


rib 




Anb/rib 


Pa/pO 


h/tf 


Ubb/tf 


Ubb/tb 


Ubf/tf 


^bf/tb 




755 
1064 


7.5 
3.5 


64 
161 


1.5 
1.5 


0.25 
0.25 


0.42 
0.74 


5 X 10"^ 
0.02 


0.39 
0.12 


5.7 
1.8 


14.6 
14.5 


5.6 
3.2 


0.77 
0.77 


23Na«Li 


1064 


5.6 


73 


1.5 


0.25 


0.93 


0.03 


0.08 


1.2 


14.7 


2.5 


0.78 



TABLE IIL Choice of optical lattice depth Vf — VJ'^ (in units of iJJ) and scattering length abf (in units of the Bohr radius 
ao) which correspond to maximal values of the supersolid transition temperature Tss, the amplitude of the supersolid density 
wave An6 (see Eqs. (1421 163|l ). and the height of the supersolid-superfluid time-of-flight peak ratio pa/po (see Eq. (|80|l ) on the 
isotropic triangular and square lattices. We consider different mixtures, optical lattice wavelengths A and bosonic fillings Ub- 
Parameters {tb,Ubb-,Uhf} (in units of tf) follow from choice of Vf,abf {Vf — 20) via Eqs. (|73l I74p . The critical superfluid to 
Mott-insulator ratio is given by Ubb/tb\sF-Mi = 26.5 (16.5) for the triangular (square) lattice. Weak-coupling analysis requires 
Xbf/tb ~ UbMoUbf /8tb < 1 for the triangular lattice. For the square lattice Mo is replaced by A^o. 



Anisotropic Square Lattice: r/ = = 0.75, nt, — 3/2 



Species 


A[nm] 


Vf 




abf 


max 


Aub/rib 


Pa/PO 


tb2/tf2 


tbl/tb2 


Ubb/tf2 


Ubb/tbi 


Ubf/tf2 ^Bf/tb 


«^Rb«K 


755 
1064 


7.5 
3.5 


6.4 
2.4 


64 
158 


0.17 
0.13 


0.34 
0.49 


3 X 10"^ 
0.01 


0.40 
0.18 


0.73 
0.52 


4.1 
1.2 


14.0 
13.3 


4.1 0.46 
2.1 0.31 


23Na«Li 


1064 


5.8 


4.7 


71 


0.17 


0.73 


0.02 


0.09 


0.62 


0.9 


15.5 


1.9 0.45 



TABLE IV: Choice of optical lattice depths Vf , Vf and scattering length abf (in units of the Bohr radius ao) which correspond to 
maximal values of the supersolid transition temperature Tss, the amplitude of the supersolid density wave Ant (see Eq. (I64|l l. 
and the height of the supersolid-superfluid time-of-flight peak ratio pa/po (see Eq. (I80|l 'l on the anisotropic square lattice. 
Other parameters follow from generalizations of Eqs. (|72l 1741) with Vf = 20, and weak-coupling analysis requires Xbf/tb ~ 
nbNoUif/8tb < 1. 



From the values given in Table HTTl and llVi we conclude 
that while it is feasible to detect the supersolid peaks 
for the square lattice geometry, they are too small to be 
detected for the triangular lattice. Therefore, another 
way to detect the density wave correlations should be 
used to confirm the supersolid nature of the system. 

This can for instance be achieved by the analysis of 
noise correlations in the absorption spectrum, where a 
density wave also leads to peaks at its characteristic 
wavevector(s) Qq, [sqI. [90|. Combined with the observa- 
tion of a (superfluid) zero momentum peak in TOF, this 
also proves the existence of the supersolid phase ^] , if 
one can exclude the coexistence of multiple phases in the 
trap. 

The coexistence of phases arises due to spatial inho- 
mogeneities in the chemical potential = /i^(;,)(x) 



introduced by an overall (harmonic) confinement. For 
example, a pure density wave, i.e. with vanishing super- 
fluid component, surrounded by a superfluid shell shows 
noise correlations similar to a supersolid, however, it does 
not show any first-order coherence peaks at Q^. As noted 
in Ref . [2^ . the differences in time-of-flight imaging be- 
tween a density wave phase coexisting with a superfluid 
shell and the supersolid phase are merely quantitative. 



VI. MOTT INSULATING PHASES 

So far we have concentrated on the case of weak- 
coupling, where the values of the interaction parame- 
ters are limited to Ubb/tb < (t^6b/4)sF-Mi, UbfMo <C 1, 
and Xbf/tb < 1, where Abf = M^U'^JUbb and tb = 
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8tb/niiUbb for triangular lattice (for the square lattice Mq 
is replaced by A^o)- The first inequality assures that the 
system is superfiuid, and not in a Mott insulating (MI) 
phase, for T > TssjTps- The second and third inequal- 
ity defines the regime where the effect of the fermions on 
the bosons can be described in second order perturbation 
theory. 

In this section we discuss the opposite region of strong 
coupling Ubb, Ubf ^ tf,tb, where the system can be de- 
scribed by an effective i-J-model. At a filling of one 
particle per site Uf + nb ~ 1, it reduces to an anisotropic 
quantum Heisenberg model. For small bosonic hopping 
if, <C t/, it turns out that the in-plane (XY) coupling is 
much weaker than the coupling of the z-components, and 
we will argue that the system has a stable, and unfrus- 
trated, antiferromagnetic ground state both on the trian- 
gular lattice for filling factors Uf =3/4, rif, = 1/4 as well 
as on the square lattice for Uf = Ub = 1/2. The fermions 
form a density wave that is characterized by the nesting 
wavevectors Qq, i.e. on the triangular lattice this is an 
antiferromagnet (AF) with a real-space Kagome-pattern, 
and on the square lattice it is the usual Neel state. For 
repulsive Ubf > 0, the bosons become localized at the 
sites where no fermion is present. This phase is exactly 
the alternating Mott insulator phase (AMI) that was de- 
scribed for the square lattice in Ref. '3^. 

We conclude that, at unit filling nf + Ub = I, a. Bose- 
Fermi mixture becomes supersolid only for sufficiently 
small interspecie interaction Ubf. It will be addressed 
elsewhere, whether the system enters a supersolid phase 
in the strong coupling regime away from unit filling, 
i.e. upon doping the AMI phase by adding or remov- 
ing bosons. Such a behavior was reported recently in one 
dimensional Bose-Fermi mixtures using quantum Monte- 
Carlo simulations |91i] . 



important property that exTp[i9pj]exp[—i9jp\ = —1. The 
purely fermionic Hamiltonian now reads 

Hb = -tb J2 [4cje'^- + h.c] -fibJ^N, 



(hi) 



Hbf ^Ubf^ miN^ , 



(82) 



where Aij = J2p^i j i^pj ~ Opi) Np. Except for the ad- 
ditional (gauge) field A^, this is the Hamiltonian of 
the two-dimensional spin- 1/2 fermionic Hubbard model, 
where //(c^) creates a spin-up (down) fermion at site i, 
and the boson- fermion interaction Ubf marks the on-site 
interaction. It is worth noting that the gauge field dis- 
appears in a one-dimensional system [935. 

In the limit of large Ubf/tf^b, it is well-known [93| that 
one can derive a t-J-model Hamiltonian that describes 
the low-energy (spin and charge) excitations of the sys- 
tem. At unit-filling Ub + Uf = 1, it reduces to the anti- 
ferromagnetic spin- 1/2 quantum Heisenberg model. 

If we follow the standard derivation (see Appendix iFl 
for details) and focus on the unit-filling case, we arrive 
at the familiar form of the quantum Heisenberg Hamil- 
tonian, except that the XY-coupling terms contain the 
( Jordan- Wigner) gauge field Aij : 



2 ^ Ubf 



2 3 



Ubf 



QZ QZ 



}~ii^f-i^b)J2s: 



(83) 



A. Derivation of quantum Heisenberg Hamiltonian 

In the limit of large Ubb, where double occupancies 
are energetically forbidden, one can replace the boson 
by spin- 1/2 operators via 6] = ^('Sj + a-iid 

Uj 1 + where s" — cr"/2, (a — x,y,z), and a" 
are the usual Pauli-matrices. We then fermionize the 
'bosonic' spins s" using the celebrated Jordan- Wigner 
transformation in two-dimensions 19211 



Cj exp 



Cj exp 



1 0pj 



iY.9„Np 



PT^j 



(81) 



1 

2 ' 



I, J, 



and —Tr<9pj<TT is the 



where Nj 

argument of the vector from site j to site p. It has the 



Here, = fjci, Sf = {rrii — Ni)/2, are proper spin 
operators, which obey [S^,S~] = 2Sf5ij. Spin-up cor- 
responds to occupation by a fermion and spin-down to 
occupation by a boson. Note that the presence of the 
gauge field reflects the different symmetry of fermions 
and hard-core bosons under exchange of two particles. 



B. Triangular lattice 

For the particular filling of n/ = 3/4, = 1/4, or a to- 
tal magnetization of {S^-) = 1/4, the ground state phase 
of the system is an alternating Mott-insulator phase for 
the bosons and a density wave phase for the fermions. 
The real-space configuration is of the Kagome-type that 
was discussed previously (see Fig.S]), only that bosons are 
now localized. Formulated in the spin-language, the sys- 
tem favors the classical unfrustrated Ising ground state 
for all values of {tb, tf}, because of the externally applied 
magnetic field in the z-direction, which is proportional to 
the number difference of bosons and fermions in the sys- 
tem. 



20 



If the system is doped away from unit filling, e.g. by 
adding or removing bosons, it can be described by an 
effective i- J- model with the additional gauge field Aij. 
It remains an open question, whether the system then 
becomes supersolid, as is the case for a one-dimensional 
Bose-Fermi mixture [9l| . 



C. Square Lattice 

Here, a stable alternating Mott-insulator phase occurs 
for double half-filling nt, = Uf — 1/2. Even without ex- 
ternally applied magnetic field (/ib — /i/), the system 
enters the (classical) Neel antiferromagnetically ordered 
ground state, because of the finite anisotropy in the spin 
coupling that occurs for > tb [H, HI, log • The XY- 
coupling is renornialized to zero, and the spins point 
along the z-axis, i.e. the system is a Mott insulator 
with a site occupation that alternates between bosons 
and fermions. This agrees with recent DMFT calcula- 
tions in Ref. [ssj . 



VII. CONCLUSIONS 

We have studied mixtures of spinless bosons and 
fermions in different two-dimensional optical lattice ge- 
ometries at fermionic fillings n f that give rise to a nested 
Fermi surface. We have shown how nesting can lead to 
supersolid formation via a density wave instability of the 
fermions. The resulting density order in the supersolid is 
characterized by the nesting vectors Q^. 

On the triangular lattice, we have thereby identified 
a novel supersolid phase with three ordering wavevec- 
tors Qi,2 = {±TT,n/V3), Q3 = (0, 27r/\/3) that give 
rise to a Kagome-pattern in real-space. We predict this 
novel phase to appear at rather low temperatures Tss ^ 
i//40 = Tf/240, and the density modulation Ani, in the 
supersolid, which is proportional to Tg^/Ubf to be weak. 
Typically, we find that the density wave only involves 
0.1% of all the bosons, which leads to Arib/nh ~ 0.05 . 

Furthermore, for temperatures such that the thermal 
energy exceeds the characteristic energy level spacing in 
the system, we have pointed out the possibility of an 
incommensurate density wave modulation in the super- 
solid phase. If thermal effects can be ignored, however, 
only the instability towards the commensurate supersolid 
remains. 

Higher transition temperatures and larger density 
wave modulations can be found for spatially anisotropic 
hopping amplitudes, or if one considers the square lat- 
tice geometry. In these cases, the nesting relation is ful- 
filled for a larger fraction of wavevectors in the first Bril- 
louin zone. We have derived transition temperatures of 
Tss — If/4 (Tf/5) for the isotropic (anisotropic) square 
lattice. The density wave now involves up to 20% of all 
the bosons. We have identified the square lattice as the 



optimal choice, since it shows the highest values of Tss 
and Arih. 

We have pointed out that introducing anisotropic hop- 
ping on the square lattice has the advantage that the 
tendency towards phase separation is weakened, which 
then only occurs above a certain inter-species interac- 
tion threshold U^f. The values of Tss and Arifc on the 
isotropic and anisotropic {rf = 0.75) square lattice are 
about the same. 

We have also predicted how to experimentally realize 
and detect the supersolid phase in two commonly used 
Bose-Fermi mixtures, *^Rb'*°K and ^'^Na^Li. The square 
lattice geometry allows for supersolid transition tempera- 
tures close to current cooling limits Tss — Tp/4 for both 
mixtures. However, since the amplitude of the density 
wave modulations grows for a smaller ratio of bosonic to 
fermionic hopping amplitudes t^/tf, we find that the de- 
tection of the supersolid phase via additional coherence 
peaks in time-of-flight absorption images becomes easier 
for smaller tf,/i/ (slower bosons). Both ^^Rb'*''K, trapped 
in a A = 1064 nm optical lattice, as well as ^'^Na^Li are 
therefore good candidates to observe supersolidity. 

Finally, we have considered the strong-coupling regime 
of Ubb,Ubf S> tf,tb and derived a quantum Heisenberg 
Hamiltonian that includes an additional gauge field due 
to the Jordan- Wigner transformation in two-dimensions. 
For filling factors of n/ = 3/4, Ub = 1/4 on the tri- 
angular and nj — Hb — 1/2 on the square lattice, the 
ground state of the strong-coupling Hamiltonian is an 
alternating Mott-insulating state (AMI) for the bosons. 
The fermions exhibit a density wave and, for repulsive 
interaction Ubf > 0, occupy all the other lattice sites. 
The order is again characterized by the nesting vectors 
Qi, which leads on the triangular lattice to the same 
Kagome-pattern we found previously for the supersolid, 
however, the bosons are now spatially localized. 
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APPENDIX A: DETAILED ANALYSIS OF THE 
LINDHARD FUNCTION ON THE TRIANGULAR 
LATTICE 

This appendix includes a detailed discussion of the 
finite temperature behavior of the Lindhard function 
x(T, q) on the triangular lattice. In particular, we point 
out that for temperatures larger than the energy level 
spacing in the system, the minimum of the function in 
k-space occurs slightly away from the nesting vectors Q^. 
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This gives rise to an instability towards a supersolid with 
an incommensurate density wave modulation. 



1. Fermion-mediated interaction 

The instability criteria are based upon analysis of the 
low temperature divergences of the fermionic Lindhard 
function x(T, q) = Y.k ^C^^ ^) "^^^^ 



FiT,q,k) 



/[e/(k),r]-/[0(k + q),T] 



C/(k)-^/(k + q) 



IT] 



(Al) 



and f{£,f,T) = [1 + exp(^//T)]-i is the Fermi function. 
The Lindhard function describes the part of the effective 
boson-boson interaction that is induced by the fermions 



U{T,ci) = Ubb + U^fX{T,q) 



(A2) 



This effective interaction is obtained after an exact inte- 
gration of the fermionic degrees of freedom using a func- 
tional integral approach, followed by a perturbative ex- 
pansion to second order in MoUbf- Therefore, this form 
of interaction is restricted to the regime of weak boson- 
fermion coupling MoUbf < 1. Here, Mq = 3/{4:TT^tf) is 
an estimate of the regular part of the fermionic density 
of states. 

At the particular filling of n/ = 3/4, the Fermi surface 
of the system both shows nesting and contains critical 
points at Q^, i.e. Vq^/(q)|q. = 0, that lead to a van 
Hove singularity in the density of states at that filling 
(see Fig. [5]). As a result, the Lindhard function diverges 
as T ^ at the wavevectors that lie on straight lines 
between q = and the three nesting vectors Qi,2.3 (see 
Fig. [T]). One can write this set of q- vectors as Af — 
{q; 3a e [0, 1] : q = aQ,}. 

At q = 0, the Lindhard function diverges logarithmi- 
cally 



X(T,0)~ -Moln 

with To = SCitf and Ci = 2e^ /tt ~ 
Euler-Mascheroni constant. 

At the nesting vectors q = Qi 
divergence is enhanced to 



To 
T 



(A3) 



1.13, where C is the 



one finds that the 



X(T,Q0 



Mo 
6 



In 



Ti 
T 



(A4) 



where Ti ^ gTq with a ~ 2.17 being determined from 
a fit to a numerical calculation of x(r, Q^) using impor- 
tance sampling Monte-Carlo integration. 

In between, for q = aQi with < a < 1, the behavior 
is different for a <C 1 and for a ~ 1, i.e. for q being close 
to the nesting vector. We investigate the two cases sep- 
arately in the following sections, and find that at finite 
temperatures, the thermal width of the Fermi functions 
in F(T, q, k) comes into play. For a finite system, how- 
ever, one can neglect thermal effects for temperatures 
below the system's characteristic energy level spacing. 



(a) 



(b) 







(d) 




FIG. 16: (Color online) Contour plot of the integrand 
F(r,q,k) < of the Lindhard function x(r, q = (1-5)Q3) = 
/ d^k F{T, q, k) for fixed temperature T = tf/W. Parts (a-d) 
correspond to different deviations of q = (1 — (5)Q3 from the 
nesting vector: 5 = 0,0.25,0.5,0.9 (a-d). Hexagon (dashed) 
denotes the first Brillouin zone. Bright regions indicate small 
absolute values of F, dark regions indicate large absolute val- 
ues. In particular, the integrand vanishes in white regions. 
For small 5 the integrand is peaked along fci — ±7r. An ad- 
ditional peak occurs at qc (Eq. (IA9p ) for nonzero S, which 
moves from Qi,2 towards Q3 along the Fermi surface. Close 
to q ~ (d), the integrand is peaked at the critical points 
±Qi,2,3 of the Fermi surface. 



2. Long- wavelength divergence 

In this section, we investigate the regime < a <C 1. 
By numerical integration, we find that the temperature 
dependence of x(2^j Q^Qi) is always logarithmically, but 
with a slope that depends on temperature. For tempera- 
tures above some a-dependent temperature it is equal 
to Mq. At the slope decreases abruptly to a slightly 
smaller value, « 2Afo/3, which holds then for T < T^. 
The transition occurs at smaller temperatures for smaller 
values of a. 

This behavior can easily be understood when the 
width of the Fermi function in the numerator of 
F(r, q, k) (Eq. (jAip ) is taken into account. For a <C 1, 
we find that the function F{T, aQi, k), if considered as a 
function of the integration variable k, is strongly peaked 
for k « Qi, because these are the critical points, which 
lead to the van Hove singularity in the density of states 
(see Fig. [HP). 

Let us consider the case of q = aQs for definiteness. 
In the limit a ^ 0, the function ^(T, aQ3,k) is equally 
strongly peaked at all six van Hove singularities. More 
mathematically speaking, if we write k = Qj + k, where 
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FIG. 17: (Color online) Shift of the contour where the de- 
nominator of F(r, q = (1 — 5)Q3,k) vanishes in k-space up 
to first order in 5 <C 1. For S = 0, the denominator vanishes 
along ki = ±7r and along the horizontal solid lines connected 
by Q3 (solid arrow). For nonzero 5, there is no shift of the 
vertical parts of the contour (to linear order in 5). However, 
the horizontal parts get shifted to larger values of k2 (hor- 
izontal dashed lines) . Whereas the numerator of F(T, q, k) 
vanishes along the lower (horizontal) part of the contour, the 
function F becomes peaked at the location qc indicated by 
the circles, since the shifted contour crosses the Fermi surface 
(connected by aQa -|- Gm (dashed arrow)). 

k is considered small, the energy denominator vanishes 
only quadratically in the small quantities (|k|,Q!). One 
finds ^/(Qi,2+k)-0(Qi,2 + k + aQ3) =i/C(a|k|), and 

0(Q3 + k) - e/(Q3 + k + aQa) - tfO{a\k\,a^) (A5) 

For finite a, the transition takes place when the 0{a^) 
term becomes of the order of the width of the Fermi 
function, which is about 2T. At this point, the peak 
of F(T, aQ3,k) at iQs is much reduced and the slope 
thus changes to 2/3 of its initial value. 

The precise terms read ^/(Qs+k)— ^/(Qs+k+aQs) = 
27r^Q!^t/ + 2y/3TTak2tf, and we can therefore estimate 

by 

TL = . (A6) 

3. Divergence at the nesting vectors 

In this section, we look at the case of a « 1. If we 
write a = 1 — 5, with small 5^1, we find that the diver- 
gent behavior of x(T, (1 — (5)Qi) changes from being of 
[In T/i/]^-type for temperatures above some a-dependent 
Ta to being proportional to InT/tf for T < Ta- Gener- 
ally, one can obtain a [InT]^ divergence of xi'^i'^) 
a wavevector K, if the Fermi surface both shows nest- 
ing (with nesting vector K) and contains critical points 
where Vq^/(q)|q;gFS = which lead to a van Hove singu- 
larity in the density of states. Additionally it is required 
though, that the nesting relation is fulfilled for the k- 
states that make up the van Hove peak in the density of 
states, i.e. that are close to the critical points. These 
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(a) (b) 




FIG. 18: (Color online) Detailed contour plot of F{T, q, k) < 
around k = Q2 for fixed q — 0.9 Q3 and various temper- 
atures logipT/t/ = -0.8,-1.1,-1.4,-1.7 (a-d). At Q2 the 
Fermi surface (light solid) has a kink and touches the bound- 
ary of the first Brillouin zone (darker diagonal line) . Brighter 
colors of the contour indicate larger absolute values of F, in 
particular, the integrand vanishes in black region (inverted 
color scheme compared to Fig. I16|l . The thermal smearing be- 
comes smaller for decreasing temperature, and the additional 
peak of F{T, q, k) is separated from the thermally broadened 
peak at the critical point Q2 for T <Tl. 

k-states have the same energy up to linear order in devi- 
ations from qj. 

Thermal smearing of the Fermi edge allows states that 
only fulfill nesting approximately, i.e. ^/(k -I- K) « 
— ^/(k), to contribute to the Lindhard integral. How- 
ever, at T = 0, it is required to fulfill the nesting relation 
exactly which demands K = Q 1,2,3 • 

Let us again specify to the case of K = (1 — 5)Q3 = 
(0, for definiteness. In Fig. [1^ we see that, for 

small 5 (parts (a-b)), the integrand F{T, K, k) is strongly 
peaked along the lines of constant ki = ±7r. The vector 
K only translates along the fc2-direction and thus pro- 
vides a mapping between two points of the Fermi sur- 
face (nesting). We also see clearly that the peaks be- 
come wider close to the critical points Qi_2i where the 
energy only varies quadratically with deviations from 
Q1.2. If we write k = Qj + k, this reads C/(Qi + k) = 

0(Q.) + O(|kp)«^/(Q,). 

To obtain a divergent behavior of type x{T, K) ^ 
— [InT/t/]^ for all T — > 0, it is required that the energy 
denominator vanishes quadratically in |k|: 

0(k)-0(k + K) = O(|k|2). (A7) 

However, this relation is only valid exactly at K = Qi,2,3- 
Slightly off, at K = (1 — (5)Qi, one finds instead that 

0(k)-0(k + K) = 5O(rk|) (A8) 
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The energy denominator now vanishes linearly with |k|. 
One can clearly observe that the peak of F{T^ K, k) close 
to k « Qi.2 narrows as 5 is increased [Fig. [TCT a-c')]. 

The plots also show the emergence of additional peaks 
close to Qi,2 [see Fig. [T6lb-c)]. which move along the 
Fermi surface towards Q3 for increasing 5. They occur 
because a part of the contour, where the energy denomi- 
nator vanishes, shifts for nonzero 5, and crosses the Fermi 
surface. This is shown in Fig. 1171 The energy denomi- 
nator vanishes along the rectangle, with vertical sides at 
ki = ±TT and horizontal sides at fc2 = n{±l + S)/V3. For 
6 = 0, the function _F(T, Q3,k) vanishes along the hor- 
izontal path, since the numerator is zero. This part of 
the path, however, shifts for nonzero d to larger ^2 values 
and crosses the Fermi surface at 



qc = (n{l-d),il + 6)7r/V3 



(A9) 



As a result, a peak in F(T, K, k) occurs around k « Qc, 
which is clearly visible in Fig. 1181 This peak is responsible 
for the fact that the minimum of the Lindhard function is 
shifted away from Qi,2,3 for intermediate temperatures. 

As long as the thermal smearing is larger than the sep- 
aration |qc — Q2I, the peak of F{T, (1 — (5)Q3, k) is actu- 
ally broader for nonzero 6 than for (5 = [see Fig. rTST a- 
b)]. However, as the temperature is lowered the sep- 
aration of the additional peak at qc from the critical 
point Q2 finally becomes larger than the thermal width 
[see Fig. [TST c-d')]. At this point, the nesting relation is 
no longer fulfilled for all the states close to Q2 (in the 
sense defined above), and the divergence of x{T,~ii) be- 
comes of type InT/tf. As a result, eventually one finds 
xC^) Q3) < x{T,K.), i.e. in the limit T — 0, the mini- 
mum of the Lindhard function occurs at Qi,2 3- 



4. Level spacing temperature Tl 

In this section we will estimate the temperature Tq, 
where the divergent behavior of xC^^ (1 ~ '5)Qi) changes 
from being [InT/tf]'^ to being InT/tf. The arguments 
we will use are similar to the ones in Sec. lA 21 

First, one can relate the thermal width of the Fermi 
function to a distance in k-space using the dispersion 
relation. Expanding the fermionic energy around the lo- 
cation of the additional peak qc (see above), yields 



C/(qc(<5) +k)~ ^6tf (ki + \/3fc2) 



(AlO) 



The thermal width of the Fermi function ±2T thus re- 
lates to k-space like 



~ 2T ~ 2T 
ki ~ ±- ^ , fc2 ~ ±- 



tfird 



y/itf-KS 



(All) 



The separation of the peak at qc from the van Hove 
singularity is given |qc - Q2I = 2^(5/ V3, we estimate 



the crossover temperature Tq to occur when 257r/\/3 
4Ta/TrStf, which leads to 



(A12) 



where 6 = 1 — a. 

For a finite lattice of Nl = unit cells, this defines a 
(level spacing) temperature Tl ~ ^/L"^, below which the 
roton gap closes at Qi,2,3, by noting that min(i5) = 2/L: 



Tr = 



tj2n^ 
L2 



(A13) 



This characteristic temperature corresponds to the spac- 
ing of energy levels in the finite system, i.e. thermal 
effects can be ignored for temperatures smaller than T^. 

It is worth noting that for the square lattice, the min- 
imum of x(T, q) always occurs at the nesting vector 
q = Qsq = (7r,7r). 



APPENDIX B: FERMION INDUCED 
INTERACTION IN REAL-SPACE 



In this section we describe the real-space form of 
the effective boson-boson interaction J7(T, q) = Uhb + 
UhfX{T, q) on the triangular and square lattice (see 
Sec. lIII A|) . We obtain the Lindhard function for wavevec- 
tors q in the first Brillouin zone by numerical Monte- 
Carlo integration. The real-space form U{T,Xi) is sim- 
ply the Fourier transform of U{T, q), evaluated at lattice 
sites X,-: 



UiT, X,) = Ubb (5x„o + U^fXiT, X,) 



(Bl) 



where (5xi,x is the Kronecker delta and the real-space 
form of the Lindhard function is given by a sum over 
wavevectors in the first Brillouin zone 



xiT, X,) 



x(T,k,)e'' 



(B2) 



The intrinsic interaction term Ubb S^i.o provides a contact 
interaction in real-space, which is repulsive for Ubb > 0. 
On the other hand, the fermion- induced part U^^xi^, x^) 
provides a long-range interaction between the bosons 
that is oscillating in sign, as can be seen in Fig. [121 

On the triangular lattice [see Fig. [TW a)] , the nearest- 
neighbor interaction is attractive, but the strongest at- 
tractive interaction is at the second nearest-neighbor sites 
±2ai_2 and ±2(ai — 3.2). The interaction is repulsive at 
all the other second nearest-neighbor sites. Clearly, this 
form of interaction gives rise to the Kagome-type density 
wave seen in Fig. |4l Note that this form of interaction 
is quite different from the usual U-V-uiodel interaction, 
which is short-ranged. There, the supersolid emerges 
for large repulsive nearest-neighbor interactions V due 
to frustration. 
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FIG. 19: (Color online) Lindhard function in real-space 
x(r, Xi) for triangular (a) and square lattice (b) at a tem- 
perature oflog^QT/tf = —2.9 (compare with Fig.|3](b)). The 
lattice site Xi = is located in the center of the lattice. The 
sign of x(r, Xi) is encoded in the shape of the lattice site: 
circles denote attractive interaction x(r, Xi) < 0, whereas 
squares denote repulsive interaction x{T,ii.i) > 0. Darker col- 
ors correspond to larger absolute values of x{T,^i)- 



On the square lattice [see Fig. [TWb)]. the nearest- 
neighbor interaction in repulsive, whereas the second 
neighbor interaction is attractive. The emergent bosonic 
density wave in the supersolid is thus of checkerboard- 
type [compare with Fig.[5]Ja)]. 

The fact that the sign of the nearest-neighbor interac- 
tion is different for triangular and square lattice reflects 
the fact that the divergence at the nesting vectors is not 
as pronounced for the triangular lattice as it is for the 
square lattice. 



APPENDIX C: CONSTRUCTION OF THE 
LANDAU-GINZBURG- WILSON FUNCTIONAL 
FOR TRIANGULAR AND SQUARE LATTICE 

Motivated by the instability analysis (see Sec. [TT|, we 
assume that the supersolid will cause condensation into 
a number of Fourier components. In addition to (bo) 
(superfluid order), we consider the possibihty of {bn.) ^ 
0, where are the wavevectors where the roton gap 
closes first. For simplicity, we focus on the case = 
Qi, i.e. on the triangular lattice we assume transition 
temperatures Tss < Tl- Here, are the nesting vectors 
of the Fermi surface. 



1. Isotropic triangular Lattice 

On the isotropic triangular lattice, there are three 
nesting vectors Qi,2,3, and the real-space wavefunction 
can be expanded as 



*(x)«^o(x)+^V.(x)e'^^' 



(CI) 



The complex fields ■0 = (V'o(x), ■0i(x), '02(x), ■!/'3(x)) are 
defined as V'i = V^L(^Qi)> where Qo = 0. The average 
is taken over a length scale much larger than the lat- 
tice spacing. They are assumed to be slowly varying in 
space-continuum and below we will eventually focus on 
the special case of completely homogeneous fields. 

The Landau-Ginzburg- Wilson free energy is con- 
structed from every possible term that is invariant under 
the symmetry group of the lattice. Under operation of 
a symmetry element of the lattice {G, u}, which acts on 
a lattice position in real space as x' — G ■ x + u, the 
real-space wavefunction ^'(x) transforms as 

^''(x') = *(G-x + u). (C2) 

In momentum space, the wavefunction transforms as 

«''(q') = *(G • q)e''i-^ " . (C3) 

The symmetry group of the triangular lattice can be 
spanned by the five generators {R, h^h, ^i, T2}, where 



R 




(C4) 



is a 3-fold rotation (choosing the origin to be a lattice 
point), and 




(C5) 



are the two inversions, and Ti,2 are the two translations 
by a Bravais lattice vector: x' = x -|- ai^2- In order to 
examine how the Fourier components ij^i transform under 
these operations, it is convenient to cast the transforma- 
tion rules into matrix form, acting on the vector space 
oi xj} — (-00, "011 "02, V'a)- This is in fact a representation 
of the symmetry group under which these momentum 
points transform. In this representation, one finds 



R 



as well as 



h 



(l o\ 
10 
1 

Vo 1 oy 



(i o\ 
0010 
0100 

yo 1/ 



(C6) 



(C7) 



i=l 



The translations take the form 

(l Q Q 
1 
0-1 

Vo 






-V 



(C8) 
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and 



T2 = 



fl \ 

0-100 
10 

\o -ly 



(C9) 



We will now shortly describe the general procedure to 
find all the quadratic and quartic terms that are invari- 
ant under the symmetry group. The N different gener- 
ators can be written as . . . , t-^"*}, (here N — 4), 

where ^Ij'' are d-dimensional matrix representations of 
the group generators, and d is the number of order pa- 
rameters considered (here d = 4). The order parame- 
ters transform under operation of a group generator like 

^'t = tifi^j and (V'O* = {tji^Y'^Pj, such that a generic 
quadratic term ^2 = J2i j ^ij'^i'^j transforms as 



^2= E A,,(t^rt^^:i.,^s,. 

ij.a.b 



(CIO) 



Invariance under all generators of the symmetry group 
requires the coefficients to obey 



A, 



ab 1 



(Cll) 



for all n — 1,...,N, which are NcP equations con- 
straining the coefficients Aij. For the quartic terms 
S'4 = J2i j k I Aijkii'ii'j4'kiph One finds that invariance 
requires 

E ^UfeKii^^^j V4c^*ld^ = ^ahcd , (C12) 

for all n = 1, . . . , iV, which are Nd^ equations constrain- 
ing the coefficients Aijki ■ 

The effective LGW free energy one derives in this way 
(see Eq. (ISGI)). reads 

J-fc^y" dx{TOo|^o|'+mi|t/>Q|2 + |aVo|'+t''|5^g|' 
2 4 

i=0 i=l 

(C13) 

where xpQ — (tAi, 7/^2, V's)- It contains ten parameters 
{mo_i, 11, uq, 1,2, 51,2,3.4} and the different terms read 

ei = |tAQr = IV'tl + IV'2|^ + IV'3l' 

+ 2 (i^inv.2p + ivinv-ap + iv-snv-ap) 

©2 - \M'\^q\' = IV^OP (IV'lP + \^2\' + 
i^l-|V^l|' + |^2|' + 1^-31' 

F2 = (V? + i4 + (V-? + ^-2 + V-a) 
-P3 = -00 ('0iV'2V'3 + cyclic) -I- c.c. 

^^4=^-0 +C.C., 



2. Anisotropic triangular lattice 

In the case of the triangular lattice with anisotropic 
hopping, one finds that its symmetry group can be 
spanned by the same generators, but without the 3-fold 
rotation R. Repeating the analysis with the reduced sym- 
metry will naturally lead to the same invariant terms as 
in the isotropic case, and some additional terms that are 
now allowed as a result of lower symmetry. The quadratic 
terms are given by [see Eq. (j65p ]: 

=mo|0o|'+™i(|0i|' + |02n+m3|V3p. (C15) 

The quartic term contains, naturally, the mass terms 
squared, the cross terms between the masses as well as 
seven additional terms. Altogether, the quartic terms 
read [see Eq. 

•^i" = ffol^ol" + + 1^-21')' + 92\H^ 

+ Ao|Vo|'(|V'l|' + IV'2|') + Ai|Vo|'|V'3|' 

+ A2|V3|'(|V'l|' + IV'2|')+Ui|0l|>2|' 

+ W2[(V'?)>2 + C.c] + U3[(V'o)*(V'? + ^2) + CC] 

+ W4[^/'SV'3^lV'2 + C.c] -I- U5[{tpl)*tpl + C.c] 

+ M6[(V'oV'l>203 + C.C.) + (1 ^ 2)] 

(C16) 



3. Square Lattice 

Only one nesting vector Qs^ = (tt, tt) occurs on the 
square lattice, so the wavefunction is expanded as ^'(x) « 
"00 -)-'0i(x)e**^=«", i.e. we consider condensation into the 
Fourier components q = 0, Qsg. The symmetry group of 
the square lattice can be generated by the translations 
along X, y, the 4-fold rotation, a reflection with respect 
to either x 01 y axis, and inversion. The ipQ component 
is invariant under all symmetry operations, and the com- 
ponent tpi only changes under the translations, where 
Tx.y : ipi = --01. 

Constructing the most general LGW free energy as 
outlined above, thus yields for the square lattice (see 
Eq. (m) 



sq,b 



d^ii 



<?o|^o|' + 5i|0i|' + 52|V'o|'|0i|'+m[0o'(V'?)*+ c.c.]}. 

(C17) 



APPENDIX D: DETAILS OF ANALYSIS OF THE 
LANDAU-GINZBURG- WILSON FREE ENERGY 
ON TRIANGULAR LATTICE 

In this appendix, we present the details of the analysis 
of the LGW free energy on the isotropic triangular lattice. 
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We begin by writing the matched parameters between the 
general LGW functional of Eq. and the microscop- 
ically derived mean- field Hamiltonian of Eq. The 
quadratic coefficients (mass terms) read 



(Dl) 



and the interaction coefficients are given by 



u 
2 



"2 



2 

2tll : 



(D2) 



as well as 



31 



52 = 



53 
54 



52 



2 = 2g 
^ 5 
2 ■ 



(D3) 



The phase diagram is obtained from minimizing the free 
energy of Eq. ([50)1 . Here, we will analytically investigate 
the free energy expresssion, and begin by identifying the 
dominant fourth order term. For this, we rewrite the 
quartic terms of Eq. ((30|) as 



-290 + 2 |V;,|2|^^f + 4F3 + 

i>]=0 



i=0 



For a fixed particle number, the first term is just a con- 
stant: ^{u — g)\tp\'^ = ^{u — g)nl, and we have to analyze 
the remaining terms only using the most general form, 



(D5) 



with Tj > 0, 00 = 0, and < < 2n for j = 1, 2, 3. 

Let us begin with g > 0. Numerically minimizing the 
last four terms of Eq. (|D4p . one finds that the global 
minimum occurs at 



ro 



nb-ri 



r2 = r-i 



0, 



(D6) 



where the last four terms sum up to {~gni,/2). Look- 
ing at the individual terms, we see that only the last 
two terms in the square brackets may favor the super- 
solid. We assume that the F3 term is dominant due to 
the prefactor of 4 multiplying it. We will find that this as- 
sumption is in accordance with the numerical result, and 
leads to the same minimum configuration. The F3 term is 
minimal when all ipi are real and have the relative phases 
(po — 0, 01,2,3 = TT, SO that 2gtpQ{ipiip2ip^+ cyclic) +C.C. — 



— 12grQrir2r^, which is minimized for ro ~ ri = r2 = 
^3 = ^/nt/'l. Writing more generally tq — ^^cosO and 
ri = r2 = = y^sin^ [0 < < it/2), we find how- 
ever, that both gJ2t>j iV'iPlV'iP = 3gr^(r^ + rf), which 
favors ri to vanish, as well as — 5|'!/'ol^i which exclusively 
favors condensation into ^/^o, hinder a supersolid configu- 
ration independently of the phases {(pi}. The remaining 
term is just a constant for tpi G R. We conclude that for 
g > 0, the superfluid has a lower free energy than any su- 
persolid phase. This is equivalent to say that a supersolid 
phase can only occur if U{T, Qi,2,3) < 0. 

We therefore turn now to the case 5 < 0. Using again 
the general form of Eq. (jD5|) . we numerically find that 
the global minimum of the last four terms of Eq. (jD4p 
occurs for a supersolid configuration 



ro = ^/mCOsOmm 

ri = r2 = rs = 



rib . n 
sm {)„ 



(D7) 



where 



1.09 rad. The system tends to order in 



a symmetric way with respect to the three nesting fields 
with the phases locked to the superfluid phase (pi.2,3 — 0, 
thus preserving the 3-fold rotational symmetry of the sys- 
tem. This is the Kagome-type order shown in Fig. 3] 
Note that O^^in denotes the lowest energy configuration 
of the fourth order term only. If we also consider the ki- 
netic (quadratic) terms, the value of the supersolid angle 
6 is generally found to be much smaller, since it costs 
kinetic energy to add a boson with wavevector Qi,2.3 to 
the system. 

If we look at the free energy term by term, and 
making the (Kagome) ansatz ro — ^/ribCOsO, ^1,2,3 — 
yjribl'ie^'^ sinfl, we find that now the last three terms 
of Eq. (|D4p favor a supersolid configuration. First, 



(D4) 25F3 = \2gror\co^(p is minimal at = and ro = ri, 
but also the two terms that were previously, for t/ > 0, 
opposing the supersolid are now favoring it as well: 
5Si>j l^iPlV'iP = 3(7r^(rg-|-r^) wants both components 
to be nonzero, and — (/IV'ol^ now opposes condensation 
into the superfluid mode and hence favors nonzero "01,2,3- 
We conclude that a Kagome supersolid phase occurs 
for a sufficiently negative g < 0. The energetically 
most favorable ordering is symmetric with respect to 
the three nesting vectors with all the phases locked: 
ipo = ^/nbcos0, ^^1,2,3 = \/ribfi sin0 with < 6* < 7r/2. 

The phase boundary between the superfluid and the 
supersolid phase can be obtained as follows. Assume that 
the system is in the superfluid phase where the fermionic 
chemical potential is given by /ij, — — 6tb -I- ribu and 
"00 — \A*b, ■01,2,3 ~ 0. Since the superfluid-supersolid 
transition is a transition between two ordered phases, we 
have to consider the influence of the fourth order terms 
onto the curvature in the ^/'i,2.3 directions by replacing 
"00 ^/rib and examine the quadratic terms in ^'1,2,3- 

We expect the transition to the supersolid phase to 
take place when the curvature in the '01,2,3-direction, be- 
comes negative. With -00 — *■ \/nb, and using the knowl- 
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edge that the minimal energy configuration is given by 
the field configuration of Eq. (|D7|) . i.e. setting = 0, 
one finds the quadratic terms m1^{rf + r'^ + r|) with 



,cff 



— —^b + 2<b + Tib [2g + u] 



(D8) 



The mass m'^ becomes negative at the critical chemical 

potential /zj;'''' = 2tb + nb[u + 2g\. Since = —QU + 

UbU in the superfluid phase, this critical value is reached 
for = "TP- The condition of changing sign 
at the superfluid-supersolid transition, then leads us to 
conclude that the bosonic system is superfluid for g > 



whereas it forms a Kagome-type supersolid for g < 



,(1) 
yc ■ 

On the other hand, if we use the fact that chemical po- 
tentials are equal at the phase boundary: fJ-^^^^ — fJ-'j^^\ 
we find that the phase transition occurs already for larger 
values of g at 



(2)^_^>9(l) 



rib 



(D9) 



Indeed, the free energy contains a third order term in 
^"1,2,3, which reads 12grorir2r3 Agy/nbrl, where we 



have replaced tq \fnb and used the Kagome-ansatz 

with = (see Eq. (|D7)) ). For values of g < gf^, this 

term leads to a local minimum at ri = ~''~'^^2c'~ ~ 
with a — 8tb + 2nbg, b = Q^/ribg and c = 4g + 3u. This 
local minimum eventually becomes the global minimum 
at (see Eq. ((55)) ) 



9c 



\2tbU 



(DIO) 



16<f, + iribu 

which marks the superfluid-supersolid phase boundary. 

APPENDIX E: DETAILS OF CALCULATIONS 
FOR ANISOTROPIC LATTICES 

1. Square lattice 

a. Density of states 

The density of states for the isotropic and anisotropic 
cases is shown in Fig. [20] and can be analytically calcu- 
lated to be 



g(M/,r/ = l)=iVoX [^l - 



(El) 



AK N + 2(F [a,k2]^F[b, fcz]) 



(E2) 



Here, k] denotes the elliptic integral of the first kind 
and K[k] denotes the complete elliptic integral of the 
first kind. The expression of Eq. (|E2p is only valid in the 



g(e) (units of No) (a) 

10 




2 4 

e (units of tf) 
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FIG. 20: (Color online) Density of states for the isotropic (a) 
and anisotropic (b) square lattice. The anisotropy parameter 
in (b) is equal to r/ = 0.75. 



region — 1 + < fif < 1 — ry, and in particular is not 
valid in the limit rj — > 1. 

It contains the quantity A'o = l/27r^i/2, which is 

a measure of the density of states away from the 

logarithmic singularity at p,f — ±(1 — r/). Also, 

for brevity we have defined the expressions feg = 

^4r//[/i2 - (rf - 1)2], k2 = y/l - k'f, where fci = 



^ — TTT, as well as tana = / Mf+''/ g^^^^ tanfe = 



Expanding the density of states around the singularity 
at jlf = — 1 + rj, we find 



g{jlf ~ 0,r/ = 1) 
9 {(j'f - {rf - l),r/ < 1) ~ Nopln 



N„\n— = Nq In (E3) 



32drf 



A/ - ('"/ - 1) 



(E4) 



for the isotropic and anisotropic lattice, respectively. 
Here, p = 1/2. /rT and -d — — , In contrast 

to the isotropic case, the anisotropic density of states is 
egular at /i/ = and equal to 



5(A/ =0,r; < 1) = 27Voif(rj) . 



(E5) 



Note that this value is significantly increased due to the 
proximity of van Hove singularities at flf = ±(1 — r/). 
This is not the case for a three-dimensional optical lat- 
tice, since the density of states is there regular every- 
where (and divergences only occur in its first derivative) . 
Therefore, supersolidity occurs at higher temperatures 
in the two-dimensional anisotropic square lattice than in 
the three-dimensional square lattice. 
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Bosonic Mean-field Analysis 



c. Fermionic Mean-field analysis 



Here, we analyze the effective bosonic mean- field 
Hamiltonian on the square lattice (Eq. ([5^ ) 



rreff 



Nl 

Q— 



6(Qa)|^a| 



1 



+ (E6) 



with 



(IV'oP + l^lP)^ 



= sin2(26l)cos2(0), (E7) 



where t/; = (^/nb cosf?, y^e**^ sin0) and Qo = 0, Qi = 
Qsq. The coefficients of the second order term read 
"T-o = 6(0) = -2(t6i + i62) - Mfc and toi = Cfc(Qi) = 
2(tbi + tvi) — fJ-b, SO from purely kinetic considerations, 
the system tends to condense solely into the superfluid 
mode -00 • However, interactions described by the pa- 
rameters (m, g) can alter the situation, and we there- 
fore minimize the zero temperature free energy density 
fb{0,(t)) = H^^ij/Nl, as a function of the angles (0,0). 
The equation 



04, fb = -gnl sin2(20) sin(20) = , 



(E8) 



can be fulfilled for three distinct cases. For 9 — 0, 7r/2 and 
arbitrary 0, which corresponds to a superfluid (SF) or a 
pure density wave (DW) phase, respectively. The third 
possibility is that = nn with integer n and arbitrary 9, 
which allows for the supersolid case ofO<6'<7r/2. At 
a minimum, it is also required that 



defb = rib sin 29 (4(^61 -I- ^62) + nbg cos 29) 

0, 



g „ 

-n^ sin cos 20 



(E9) 



which can again be fulfilled by 61 = 0, 7r/2 (SF, DW) for 

-2(tbl-l 
2nhg 



arbitrary 0, or = rnr and cos^ 9ss - ""^"^''"^^'"^ 



which corresponds to a supersolid, if < cos^ 9 < 1. 

If we compare the energy of the three cases, we im- 
mediately find that the DW always has larger energy 
than the superfluid fb{9 = 7r/2,0) — fb{9 — 0,0) — 
4rtb(i6i +tb2) > 0. However, if we compare the energy of 
the superfluid with the supersolid, we find that the global 
minimum can occur at 9ss only for sufficiently negative 
g <0. Specifically, 



/,(O,0)-/,(0ss,«^) = 



{nbg + 2{tbi+tb2))^ 
25 



(ElO) 



is positive for negative g. However, in order for the analy- 
sis to be self-consistent, it is required that 1 > cos^ 9ss > 
0, which for 5 < requires that g < —2{tbi -\- tb2)/nb. 
The superfluid to supersolid phase boundary, which is 
defined by fb{0,(j)) — fb{&ss,n.TT), occurs at the critical 
interaction strength 



9sqA 



2(tf)l + tb2) 



Ub 



(Ell) 



with the supersolid occuring for g < (?sq,c- 



We start from the Hamiltonian H = ijj^^ ^ ^ 
similar to Eq. ([55]) . where the bosonic operators bo Q^ are 
replaced by the complex fields "00, 



H 



(1) 



= 4(4i + tb2)\A\^ + ^[nl + (^o0t + c.c.) 



Nl ' ' 2 



k a,P 



C/(k) Ubfii^ori+c.c.) 



UbfiipQipl + C.C. 



-?/(k) 



(E12) 



where the nesting relation ^/(k+Qi) — — ^/(k) was used, 
and the sum over wavevectors is restricted to 1/2 of the 
first Brillouin zone. With = ro, 4'i = riexp(z0) and 
defining A = 2Ubfrori, the fermionic eigenenergies read 



S(k,A)±=±./e/(k)2 + A2cos2 



(E13) 



We identify A as the emerging gap in the fermionic spec- 
trum. 



2. Triangular lattice 



Here, we analyze the effective Hamiltonian of Eq. ([65 



Nr 



+ -{u + vV + gW)\ip\'^ ., 



(E14) 



where m, = ^b{Q^), u ^ C/(0), v = C/(Qi,2), g = f/(Q3) 
and 

viv^r = 2 + m^) + i^3i'(i0ii' + 1021')] 

-I- [0o01'02V'3 + '0oV'l'02'03 + 2000^-02 "03 + C.C.] 
= (0003*+ 1^102* +C.C.)2. 

(E15) 

Since < V < 2 and < VV < 1, stability requires 
u + 2v> \g\ for u, w > 0, g < 0, where we anticipate that 
the supersolid only occurs for negative g. The Lindhard 
function is regular at Q1.2, such that t; > u. In particular, 
V is positive, since with x(2^5 Q1.2) ~ —Mq, one finds 
v = Ubb + U^fXiT, Qi,2) > 0, because MoU^f/Ubb < 1. 

The masses read explicitly mo = — /i^ — 2{tbi + tb2), 
TOi = — /ife -I- 2tbi and = — /i^ — 2(tf,i — 2tb2), so from 
purely kinetic energy considerations, the system prefer- 
ably condenses into the 0o-mode only, i.e. is superfluid. 

With this is mind, we turn to analyze the interaction 
terms: as w > it is energetically favorable to minimize 
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V, which is achieved by either having the pair of fields 
(■007 V's) vanish or the pair (^/'i,V'2)- For positive 5 > 0, 
also the W-term is minimized for either ■00 = "03 = or 
■01 = ^2 = 0. As a result, kinetic energy considerations 
will select the superfluid state for positive g. 

Let us turn to the case of negative g. We observe 
that in the W-term, the field V'3 couples to the superfluid 
component "00. The system can therefore, possibly, lower 
its energy, compared to the superfluid, by allowing for 
nonzero ^3 while still having ■0i = ^02 = such that 
V = 0. In this subspace of possible field values, the 
mean-field Hamiltonian is of the same form as on the 
square lattice (see Sec. IIV 



APPENDIX F: DETAILS OF DERIVATION OF 
QUANTUM HEISENBERG HAMILTONIAN 

We start from the purely fermionic Hamiltonian of 

Eq. im) 



(i.3) 



Hbf ^Ubf^ m^N^ , 



(Fl) 



where Ai. 



E 



ditional (gauge) field Aij, 



Except for the ad- 
this is the Hamiltonian of 
the two-dimensional spin-1/2 fermionic Hubbard model, 
where /|(c|) creates a spin-up (down) fermion at site i, 
and the boson- fermion interaction Ubf marks the on-site 
interaction. 

In the case of unit filling, the ground states of the ze- 
roth order (in tf_b/Ubf) Hamiltonian U = UbfJ2i''^i^i 
are states where each site is occupied by exactly one par- 
ticle, either a fermion or a boson. Since all other states 
involve at least one doubly (or multiply) occupied site, 
one can divide the Fock space into T = S ^V, where S 
is the degenerate ground state manifold of U, and T) con- 
tains all other states. Defining the projection operators 
onto the two subspaces Ps and , one can partition the 
Hamiltonian H = Hm — Eil/^fc-^i + ^^ffni) into 



^^^^{Ps{T + U)Ps PsTPd 



(F2) 



y PdTPs Pd{T + U)Pd) 
where T^Tb + Tf = - E(.,,> Mc,e'^'^ +h.c.) 



contains the hopping terms, that we will treat in pertur- 
bation theory, and Hm denotes the Hamiltonian for a 
fixed particle number difference (or magnetization) . The 
effective Hamiltonian H"}^, acting in the ground state 
manifold S only, can be obtained by 



Using (^ B)" 
Hamiltonian 



Ps{E - H)-^Ps = [E- Hfl{E)]-^ 



(F3) 



[A-BD-^C)-^, one finds the effective 



H'i^ ^ PsTPs + PsTiPoiE - (T + U)]PDr^TPs 



PsTPl 
Ubf 



m^N^PoTPs + 0{tljlUlf,E/Ulj) , 

(F4) 



where we have used that PsTPs ~ for unit filling, and 
we have expanded to lowest non-trivial order in tffi/Ubf- 
For T = Tb+Tf, we obtain four terms H^f = a'+ B + 
C + D, which read 



A 



B 



C 



tbtf 
tbtf 



cpe 



''.Aal3 ft f 

J gja 



D = 



' TT E/ fafpflfa' 



(F5) 
(F6) 
(F7) 

(F8) 



Defining proper spin operators via 5*+ 
(toq, — Na)/2, we arrive at Eq. 



fl' 



a 1 



Sp e 



iAb 



S^Spe 



ntl+tj)r 



u, 



bf 



QZ QZ 



} - (m/ - Mb) E 



(F9) 



Here, spin-up corresponds to occupation by a fermion 
and spin-down to occupation by a boson. This Hamil- 
tonian is of the familiar form of the spin-1/2 quantum 
Heisenberg Hamiltonian, however, it contains the addi- 
tional Jordan- Wigner gauge field Aij reflecting the differ- 
ent symmetry of hard-core bosons and spinless fermions. 
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